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Preliminary  Analysis 


1.1  Introduction 

The  earth’s  atmosphere  can  affect  the  propagation  of  laser  beams  significantly.  Large-scale  struc¬ 
ture  alters  the  beam  trajectory,  while  the  cumulative  effect  of  small-scale  structure  induced  by 
atmospheric  turbulence  imparts  a  random  modulation  that  can  defocus  the  beam.  These  phenom¬ 
ena  are  accommodated  in  propagation  models  by  partitioning  the  refractive  index  field  as 


n(r,  t)  =  n(r,  t)  +  Sn(r,  t), 


(1.1) 


where  n{r,t)  represents  the  background  and  Sn{r,t)  represents  small-scale  structure.  The  parti¬ 
tioning  is  somewhat  arbitrary,  but  Jn(r,  t)  can  be  constructed  formally  by  high-pass-filtering  n(r,  t) 
on  the  spatial  variable.  The  performance  of  the  optical  system  in  the  Airborne  Laser  Extended 
Atmospheric  Characterization  Experiment  (ABLE  ACE)  is  not  affected  by  structure  larger  than 
a  few  kilometers;  moreover,  the  temperature  probe  that  measures  the  turbulent  structure  scans 
rapidly  enough  that  variations  within  typical  analysis  segments  can  be  ignored.  Thus,  high-pass 
filtering  the  output  of  the  temperature  probe  prior  to  digitization  effectively  extracts  the  fluctuat¬ 
ing  component.  The  filtered  temperature  fluctuations  determine  the  refractive  index  fluctuation, 
Sn{'Vt,t),  via  the  relation 


Sn  —  — 


79  X  10-®p 

r2 


ST. 


(1.2) 


Temperature,  formally  a  passive  scalar,  fluctuates  in  response  to  the  velocity  turbulence.  If  the 
turbulence  has  achieved  a  sustained  homogeneous  and  isotropic  state,  the  power  spectral  density 
(PSD)  of  the  temperature  fluctuations  admits  an  inertial  subrange  with  a  5/3  power-law  spectral 
index.  In  a  real  environment,  however,  the  probe  encounters  turbulent  structure  in  various  states 
of  evolution  and  decay.  Thus,  the  ideal  model  applies  only  to  disjoint  data  segments,  potentially 
of  varying  length.  As  a  consequence,  any  attempt  to  measure  parameters  that  characterize  the 
turbulence  will  fluctuate  in  response  to  sampling  errors  as  well  as  in  response  to  processes  intrinsic 
to  the  structure  evolution.  Careful  consideration  must  be  given  to  data  segmentation,  parameter 
estimation,  and  the  analysis  of  the  parameter  fluctuations  to  properly  reconstruct  the  turbulence 
for  modeling  purposes. 

This  report  demonstrates  new  analysis  procedures  that  treat  the  standard  power-law  parameters 
as  random  processes  that  have  a  richer  structure  than  mean  values  corrupted  by  sample  noise. 
Ultimately,  the  intrinsic  structure  of  the  parameter  processes  may  provide  new  insight  into  the 
physics  of  the  evolving  turbulence.  In  the  near  term,  the  analysis  is  important  for  determining 
the  structure  and  deployment  of  phase  screens  used  to  emulate  the  effects  of  the  refractive  index 
structure  in  simulations  of  the  ABLE  ACE  optical  processor. 

Spectral  analysis  has  played  a  seminal  role  in  the  theory  and  measurement  of  turbulence,  al¬ 
though  the  guidelines  for  standard  spectral  analyses  methods  are  derived  from  statistically  homoge¬ 
neous  or  locally  homogeneous  models.  Within  the  context  of  these  procedures,  a  fixed  sample  rate 
determines  the  smallest  spatial  structure  that  can  be  resolved  (Aa;  =  V At),  but  there  is  consider¬ 
able  flexibility  in  choosing  the  data  segmentation  and  subsequent  parameter  extraction  operations. 
In  particular,  data  segmentation  determines  the  minimum  scale  at  which  the  parameters  that  char¬ 
acterize  the  turbulence  can  vary.  Thus,  it  is  important  to  develop  procedures  that  cleanly  separate 
intrinsic  structure  in  the  extracted  parameters  from  sample  noise.  The  research  described  in  this 
report  addresses  this  problem  within  the  context  of  the  ABLE  ACE  objectives. 

In  Part  I  we  describe  the  procediues  that  have  been  developed  to  isolate  a  power-law  subrange 
in  the  temperature  fluctuation  data.  These  procedures  are  reviewed  within  the  context  of  standard 


1 


parameter  estimation  procedures  that  have  been  applied  to  the  ABLE  ACE  data.  This  is  important 
because  the  power-law  range  of  spectral  wave  numbers  that  represent  scale-invariant  turbulence 
varies  with  the  strength  of  turbulence.  While  this  suggests  using  a  variable  window,  it  is  desirable 
to  fix  the  wavelength  range  over  which  turbulence  parameters  are  extracted.  The  low-frequency 
limit  is  determined  by  the  block  length  for  analysis  as  described  below.  The  upper  limit  ideally 
would  be  determined  by  dijBFusion  processes.  Vagaries  of  the  data  further  restrict  the  spectral  regime 
where  viable  parameter  estimates  can  be  made.  The  analysis  presented  in  Part  I  highlights  the 
variability  of  the  the  power-law  slope  and  strength  parameters  that  ultimately  determine  system 
performance.  In  Peirt  2  these  parameter  variations  are  subjected  to  an  in-depth  statistical  analysis 
that  provides  a  new  framework  for  interpreting  the  turbulence  parameters. 

1,2  Data  Processing 
1.2.1  Data  Extraction 

The  preliminary  processing  of  the  aero-thermal  probe  data  follows  the  procedures  specified  in  the 
ABLE  ACE  report  PL-TR-96-1084,  Part  1,  Chapter  6.  The  raw  data  have  been  calibrated  to 
represent  temperature  fluctuation  in  degrees  Celsius.  Although  the  specifics  of  the  mission  (date, 
location,  altitude,  airspeed,  etc.)  are  not  contained  in  all  the  data  files  made  available  to  us,  the 
approximate  values  summarized  in  Table  1.1  have  been  used  to  convert  temperature  fluctuations 
to  refraction  index  fluctuations  via  (1.2)  when  the  actual  values  were  not  known.  The  Nyquist 
frequency  corresponding  to  the  2-cm  sample  separation  is  —  w/Ax  =  157  rad/m.  Following 
the  ABLE  ACE  terminology,  a  block  of  length  Lb  containing  Nb  =  Lb /Ax  samples  is  the  smallest 
segment  to  which  spectral  analysis  is  applied.  Block  sizes  are  constrained  so  that  Nb  =  2^.  Data 
segments  are  generated  by  averaging  the  results  from  several  blocks.  The  block  length  determines 
the  finest  level  of  segmentation  and,  thereby,  the  largest  spatial  scale  that  can  be  resolved,  while 
the  segment  length  determines  the  smallest  scale  over  which  parameter  estimates  can  vary. 

In  building  data  segments,  overlap  (blocks  common  to  more  than  one  segment)  is  avoided  be¬ 
cause  it  introduces  correlation  in  the  fluctuation  of  the  estimated  parameters  that  complicates  their 
interpretation,  as  described  in  Part  2.  For  the  same  reason,  it  is  desirable  to  estimate  parameters 
prior  to  any  averaging  operations  that  are  not  intrinsic  to  the  spectral  estimation  procedures.  Be¬ 
cause  of  the  fluctuation  statistics  of  standard  PSD  estimates  and  the  nonlinear  logarithm  operation 
that  is  applied  prior  to  least-squares  extraction  of  the  turbulence  parameters,  however,  some  com¬ 
promises  are  necessary.  We  show,  nonetheles,  that  wavelet  scale  spectra  are  less  prone  to  such 
bias. 


Table  1.1:  Parameter  Values 


Parameter 

Value 

Comment 

Altitude 

Temperature 

Pressure 

40,000  ft 
210°  K 
200  mb 

Prom  Table  6-1 

From  Table  6-1 

From  Standard  Atmosphere 

Sample  separation 

2  cm 

*^=12  kHz  @  210  m/s 

1.2.2  Parameter  estimation 

Turbulence  is  characterized  mainly  by  the  inertial  subrange  where  the  spectral  density  function 
varies  as  CK~^.  Our  objective  is  to  obtain  stable  unbiased  estimates  of  C  and  rj,  which  we  interpret 
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as  random  processes  in  their  own  right.  A  standard  approach  applies  a  linear  least-squares  fit  to 
where  ^(Kn)  is  an  estimate  of  the  PSD  at  spatial  wave  number  nAK.  There  are  a 
number  of  potential  problems  with  this  approach.  The  location  of  the  inertial  subrange  is  unknown 
a  priori;  moreover,  the  logarithm  of  the  spectral  estimate  has  complicated  fluctuation  statistics  that 
can  introduce  biases  in  their  own  right.  To  provide  an  objective  criterion  for  locating  power-law 
segments,  the  least-squares  fit  has  been  determined  over  spectral  windows  starting  at  the  lowest 
non-zero  frequency  and  continuing  to  successively  higher  frequencies  until  the  Nyquist  frequency  is 
encountered  at  the  upper  end  of  the  fitting  interval.  A  uniform  power-law  segment  should  yield,  on 
average,  the  same  fit  parameters  irrespective  of  the  fit  range.  Beyond  that,  we  rely  on  simulations 
to  verify  that  the  sampling  noise  introduced  by  the  parameter  extraction  process  is  being  isolated 
and  removed. 

Standard  processing  starts  with  an  estimate  of  the  spectral  density  function  (PSD).  The  un¬ 
normalized  discrete  Fourier  transform  (DFT)  evaluated  by  most  FFT  algorithms  is  defined  as 

Nb-I 

Xn  =  Xf.exp{~27rink/NB}^  (1.3) 

The  following  relation  between  the  sample  second  moment  of  the  process  and  the  total  spectral 
intensity  can  be  verified  by  direct  computation: 

Nb-I  Nb~1  ,  . 

wT,  /r^B,  (1.4) 

k=zO  n=0 


Thus,  the  “raw”  periodogram, 


PSDn  =  \Xn  /Nl 


(1.5) 


is  the  most  basic  PSD  estimator.  For  the  class  of  real  processes  of  interest  here,  the  PSD  is 
symmetric.  It  follows  that  the  Parseval  relation  that  equates  the  sample  second  moment  to  the 
total  PSD  intensity  can  be  written  as 


Nb-I 
^  k-Q 


xl 


JVs/2-l 

PSDq  +  2  ^  PSDn  +  PSDf^^/.^. 

n—\ 


(1.6) 


It  is  well  known  that  the  raw  periodogram  has  fluctuation  statistics  that  do  not  converge  with 
increasing  block  size.  To  circumvent  this  problem,  one  uses  a  combination  of  weighting  functions 
and  block  averages.  It  is  sufficient  for  our  immediate  purposes  simply  to  note  that  the  weighting 
function  is  normalized  to  minimize  the  bias  introduced  by  the  spectral  smoothing  that  results  from 
applying  the  weighting  function.  Other  operations  such  as  mean  removal,  detrending,  and  zero 
padding  further  complicate  the  interpertation  of  the  parameter  fluctuations.  Thus,  we  have  not 
employed  any  techniques  beyond  applying  standard  window  functions  and  block  averaging. 

For  ABLE  ACE  it  is  also  important  to  maintain  absolute  scales.  Upon  extracting  a  parametric 
representation  of  the  form  ^{K^)  ~  CK~'^^  it  is  desirable  to  present  C  in  units  such  that  a 
statistically  equivalent  realization  of  the  process  will  reproduce  the  structure  sensed  by  a  light 
beam  traversing  the  same  path.  Analogous  to  (1.6),  the  defining  relation  for  the  one-dimensional 
probe  data  is 

=  (1.7) 

where  the  angle  brackets  denote  an  ensemble  average.  That  is,  a  realization  of  the  process  based  on 
estimates  of  its  underlying  spectral  characteristics  should  preserve  the  relation  between  the  second 
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moment  and  the  integrated  spectral  intensity.  For  uniform  spatial-frequency  sampling,  the  proper 
scaling  is  27rPSDn/AK.  In  Part  1  of  PL-TR-96-1084,  Equation  (6-11)  defines  (after  an  erroneous 
square  root  has  been  removed)  the  one-sided  temporal-frequency  PSD 

_  82|rfZ(/)|2 

®  -  3—Aj—’  ‘  ’ 

where  A/  =  1/{N At)  and  the  factor  8/3  provides  the  renormalization  that  preserves  (1.6)  for  the 
Hanning  window.  In  our  notation,  \dZ{fn)\‘^  ^  PSDn^  The  frequency-interval  scaling  defines  the 
units;  however,  there  is  no  convention  for  indicating  radian  units  for  spatial  frequencies.  Radians 
per  meter  will  be  used  throughout  this  report. 


1.2.3  Scale  spectra  versus  conventional  spectra 

There  is  considerable  interest  in  wavelet  scale  spectra  as  an  alternative  to  the  DFT-based  estimation 
just  described.  The  data-probing  power  of  wavelet-based  analysis  is  discussed  in  Part  2.  Here,  we 
concentrate  on  the  relation  between  the  scale  spectrum  derived  by  averaging  the  wavelet  detail 
coefficients  over  the  delay  variable  and  the  periodogram  estimate.  This  type  of  analysis  is  usually 
developed  by  usingthe  scaled  “mother”  wavelet.  Here  we  present  a  complementary  analysis  based 
on  the  discrete  wavelet  transform  as  derived  from  the  low-pass  and  high-pass  filter  coefficients  that 
generate,  respectively,  the  wavelet  scale  function  and  the  wavelet  function  itself. 

To  set  up  the  problem,  let  represent  a  sequence  of  discrete  random  variables.  From  Section 
1.2.2,  the  raw  periodogram  estimate  of  the  PSD  is 
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x{n)f 


where 

x{n) 

is  the  DFT  of  x^.  For  a  homogeneous 


N-l 

—  ^  exp{— 27rwA;/Ar} 

k=Q 

process,  x^  admits  the  spectral  representation 


pKn 

Xk=  exp{iKkAx}  dZ{K)^ 
J~Kn 


(1.9) 


(1.10) 


(1.11) 


where  dZ{K)  has  the  formal  orthogonal-increments  property 


{dZ{K)dZ%K^))  -  6{K  -  K')^K)/{27r).  (1.12) 


Without  loss  of  generality,  we  can  assume  either  that  the  spectrum  has  no  content  beyond  Kn  ^ 
n/Ax^  or  that  any  content  beyond  the  Nyquist  range  has  been  “aliased”  into  the  Nyquist  range. 
Using  the  orthogonal-increments  property  followed  by  some  straightforward  manipulations,  it  can 
be  shown  that 


sin^  [t(^  ~  nAK)Ax^ 
sin^  [i(K  -  nAK)Ax\ 


$(K) 


dK 

2-17 


^nAK), 


(1.13) 


where  AKAx  =  277 fN.  The  approximation  holds  with  increasing  N  for  n  fixed  at  any  value  within 
its  range.  If  a  properly  normalized  weighting  function  had  been  used,  a  similar  relation  involving 
a  more  complicated  spectral  smoothing  operation  would  result.  Either  way,  for  the  purposes  of 
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calculating  averages  of  quantities  involving  periodogram  estimates,  the  DFT  coefficients  admit  the 
following  approximate  orthogonal-increments  property: 

{x{n)x*{n'))^N{Pn)Sn-n',  (1.14) 

Effectively,  we  are  assuming  that  any  bias  introduced  by  the  weighting  function  is  negligible  for 
the  class  of  power-law  spectra  that  are  being  investigated. 

We  now  consider  the  multi-scale-filtering  representation  of  the  discrete  wavelet  decomposition. 
If  Si  and  hi  represent,  respectively,  the  average  and  detail  filter  coefficients  associated  with  a 
wavelet  scale  function,  the  following  recursion  defines  and  The  latter  is  the  discrete  wavelet 
decomposition  at  scale  j: 


~  for  fe  =  0,  •  *  * ,  iV  —  1  where  N  =  2^. 
For  j  =  1,  •  *  • ,  M  with  A;  =  0,  •  •  • ,  2^^^  —  L  H- 1 


ai  = 


4 


L-1 

1=0 

L-1 

(1.15) 

52 

1=0 

(1.16) 

From  the  inverse  of  (1.10),  it  follows  that 


I 

4  =  H  exp{27rmfc/iV}.  (1.17) 

n=0 

Substituting  (1.17)  into  (1.15)  and  (1.16)  for  j  =  1  yields 

I  =  ^  I  |«(n)exp{27rmA:/(iV/2)},  (1.18) 

where 

=  ^  I  I  exp{— 27rOT//iV}exp{27rm(L  —  1)/Ar}.  (1-19) 

To  obtain  a  relation  that  contains  only  discrete  spectral  domain  quantities,  we  substitute  (1.18) 
into  (1.15)  and  (1.16).  For  j  =  2,  the  “new”  spectral  representation  of  the  weighting  coeflBcients 
effectively  involves  an  N/2  point  DFT  extended  to  N  by  periodic  repetition.  Thus,  by  induction, 
for  i  =  2, .  • .  ,M  with  fe  =  0,  •  •  •  ,2^“^'  -  {L  -  1)  +  {L  -  2)/2^ 

.  1^-1^ 

4  =  ITF  Z)  S{n]j)^n)exp{2wink/2^  (1.20) 

n=0 


where 


i-2 

S{n-J)  =  d{n/2^~'^)  ]][  ^n/2‘). 


l=zQ 


(1.21) 


Equation  (1.20)  describes  the  recursion  leading  to  the  detail  coefiicients  as  a  scale-dependent  filter¬ 
ing  operation.  In  (1.21),  the  product  term  represents  the  successive  low-pass-filtering  and  down- 
sampling  followed  by  repetition  of  the  surviving  coefficients  to  maintain  N  terms. 
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The  wavelet  scale  spectra  are  obtained  by  forming  the  sample  average 


1  ■ 

^  k=^l 

where  Nj  represents  the  number  of  detail  coeffecients  available  at  scale  j.  The  result  is  indpendent 
of  fc,  but  difficult  to  interpret.  To  the  extent  that  (1.14)  is  valid,  however,  it  follows  that 


'  •  71=0 

In  the  same  way  that  the  standard  spectral  weighting  is  used  to  extract  an  estimate  of  the  PSD,  if 
is  concentrated  around  2  the  scale  spectrum  should  be  proportional  to  the  power 

spectral  density  evaluated  at  2'~'^^^AK.  Figure  1.1  shows  the  filter  functions  through  scale  7  for  the 
Haar  wavelet.  The  first  four  scales  are  superimposed  and  plotted  against  a  uniform  frequency  scale. 
The  remaining  wavelets  are  plotted  individually  against  a  logarithmic  frequency  scale  to  show  the 
low-frequency  behavior  more  clearly.  Each  filter  has  its  peak  at  although  the  j  =  1  filter 

is  half-band  in  its  frequency  response.  For  comparison,  Figure  1.2  shows  the  corresponding  filters 
for  the  db5  wavelet. 


Figure  1.1:  Spectral  filter  functions  for  Haar  wavelet  with  1024  points.  Only  the  positive  frequencies 
are  shown,  and  each  filter  is  normalized  to  unity  at  its  peak  value. 

Because  of  the  complexity  of  the  scale  filters,  it  is  difficult  to  make  a  general  statement  about 
the  precise  relation  between  the  expectation  of  the  scale  spectrum  and  the  expectation  of  the 
periodogram.  To  investigate  this  relation,  we  have  made  comparisons  between  PSD  and  wavelet- 
scale  spectra  applied  to  simulated  data.  We  have  also  evaluated  the  wavelet  spectral  filter  functions 
for  the  theoretical  form  of  the  spectral-density  function  used  to  generate  the  simulated  data  to  which 
the  wavelet  and  PSD  analyses  were  performed.  Figures  1.3  and  1.4  show  wavelet  scale  spectra  results 
for  the  dbl  (Haar)  and  db5  wavelets.  With  Haar,  the  measured  wavelet,  the  theoretical  average 
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Spectral  Weighting  for  dbS  wavelet  (scalespcRevI  m) 


Figure  1.2:  Spectral  filter  functions  for  dbS  wavelet  for  compeirison  to  Haar  wavelet  shown  in  Figure 
1.1. 

scale  spectra  estimate,  and  the  theoretical  spectra  agree  quite  well.  For  the  db5  wavelet  there  is 
a  scale-dependent  bias  between  the  theoretical  and  measured  scale  spectra  that  is  not  reproduced 
by  the  theoretical  computations.  If  the  theoretical  computation  of  the  filter  coefficients  is  in  error, 
it  is  possibly  due  to  an  inadequate  representation  of  the  filter  functions  at  the  logarithmic  sample 
intervals.  The  filter  functions  for  the  wavelet  at  scale  nine  is  superimposed  on  each  plot.  We  see 
that  the  sidelobes  are  not  critically  sampled. 

To  complete  our  evaluation  of  the  parameter  estimation  procedures,  we  compare  the  parameter 
estimates  as  extracted  from  log-linear-least-squares  fits  applied  to  standard  PSD  estimates  and  to 
Haar  wavelet  scale  spectra.  In  all  cases,  the  biases  in  the  slopes  were  found  to  be  quite  small.  The 
strength  parameter,  however,  is  more  complicated.  The  PSD  estimates  exhibit  a  bias  that  varies 
systematically  with  the  number  of  blocks  averaged  prior  to  applying  the  parameter  extraction.  The 
overall  bias  also  depends  on  the  weighting  function.  This  characteristic  is  evidently  known,  and 
it  is  clear  that  it  is  caused  by  the  complicated  statistical  fluctuations  introduced  by  the  nonlinear 
logarithmic  operation  applied  to  the  spectral  estimates  prior  to  the  linear  least-squares  curve  fitting. 
Figure  1.5  summarizes  the  results  for  both  the  PSD-based  and  Haa.r  wavelet  scale  spectra  results. 
The  Haar  wavelet  scale  analysis  produced  a  small  (less  than  0.2  dB)  fairly  uniform  bias  up  to  eight- 
block  averages.  Beyond  that  the  available  range  of  unambiguously  sampled  wavelengths  evidently 
begins  to  affect  both  the  standard  and  the  wavelet  scale  spectra.  Because  of  the  desire  to  maintain 
high  resolution  with  as  little  bias  as  possible,  the  Haar  wavelet  scale  spectra  were  used  for  the 
analysis  of  the  parameter  fluctuation  process  as  presented  in  Part  2. 
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Circular  frequency,  K  (radiarts/m) 


Figure  1.3:  Comparison  of  periodogram  and  scale-spectrum  PSD  estimates  for  dbl  wavelet 


Figure  1.4:  Comparison  of  periodogram  and  scale-spectrum  PSD  estimates  for  db5  wavelet 
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Bias  effects  -  averages  over  1 0  simulated  realizations  (readnplotavg.m) 


Figure  1.5:  Comparisons  of  strength  parameter  estimates  from  standard  spectral  estimators  and 
wavelet  scale  spectra. 
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1.3  Analysis  of  file  95135-04  data 


Figure  1.6:  Overview  of  temperature  fluctuation  data  from  file  95135-04. 

The  data  from  file  95135-04  have  been  targeted  for  in-depth  analyses  by  different  research 
groups.  Figure  1.6  shows  the  data  set  in  its  entirety  but  decimated  by  25  samples  for  more  eflicient 
display.  The  data  show  structure  on  varying  scales  up  to  the  4-km  large-scale  (low-frequency) 
cutoff  of  the  temperature-probe  filter.  Irrespective  of  any  other  consideration,  we  observe  the 
abrupt  increase  in  the  intensity  of  the  small-scale  structure  at  distances  beyond  60  km.  To  capture 
the  transition,  a  segmentation  length  of  one  kilometer  is  needed.  To  probe  the  structure  at  low 
spatial  frequencies.  Figure  1.7  shows  a  PSD  obtained  from  the  entire  decimated  data  set  shown 
in  Figure  1.6.  A  prominent  feature  is  the  low-frequency  enhancement  beginning  near  the  low- 
frequency  filter  cutoff.  The  enhancement  extends  downward  to  spatial  frequencies  approaching  27r 
radians/km.  The  wavelet  analysis  presented  in  Part  2  isolates  this  structure,  which  we  believe  to 
be  spurious,  but  for  our  piu-poses  here  we  simply  note  that  the  spectral  components  approaching 
27r  X  10“^  radians/meter  are  suspect  and  should  be  avoided  for  parameter  estimation. 

We  also  note  the  narrow  spectral  line  at  the  spatial  frequency  corresponding  to  1.2  m.  This 
feature  is  also  suspect,  but  we  shall  see  that  it  is  well  outside  the  range  where  extraction  of 
turbulence-related  parameters  is  meaningful. 

1.3.1  Analysis  by  data  sub- segment 

Based  on  the  preliminary  survey  of  the  data,  PSD  estimates  were  performed  with  block  sizes  of 
2^3  =  (163.8  m),  2^4  (327.6  m),  2^^  (655.36  m)  and  2^6  (1.31  ^m).  Figure  1.8  shows  the  sliding 
least-squares  estimates  of  the  power-law  index  and  structure-constant  parameters  obtained  from 
the  first  50  km  of  one  10-km  segment  within  the  data  set.  Each  data  point  corresponds  to  the 
one-decade  estimate  centered  on  the  spatial  wave  number  at  which  it  is  plotted.  The  data  fits 
are  obtained  from  averages  of  the  sub-segment  PSD  estimates.  The  longer  block  lengths  provide 
more  low-frequency  information  but  poorer  stability,  since  fewer  blocks  are  available  for  averaging. 
The  horizontal  dashed  line  in  the  upper  frame  corresponds  to  the  Kolmogorov  value.  In  the  lower 
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ABLE  ACE  DecimateBy20,  nblock  =  65536,  psd  avg  over  4  block  s  (CompReport.m) 


I - .  .  ...»l - ■  I  ■  mill _ ■  . 1  . 1  •  ■  . . I  . I  . I 

10"®  10"^  10"^  10"^  10"'  lo'’  10'  10^  10^ 

Circular  Spatial  frequency  K  (rad/m) 


Figure  1.7:  PSD  estimate  derived  from  20:1  decimation  of  file  95135-04  data  after  conversion  to 
refractive  index  fluctuation  units. 

frame,  the  dashed  line  indicates  the  “clear”  level  at  which  significant  system  effects  can  occur. 
The  spectra  all  show  shallowing  (less  negative  index  parameters)  at  /C  =  2  rad/m  followed  by  a 
steepening  at  higher  spatial  wave  numbers.  The  dependence  of  the  results  on  segmentation  shows 
that  the  power-law  characteristics  in  the  low-frequency  range  are  not  well  defined. 

The  results  of  applying  the  same  parameter  extraction  procedure  to  the  second  50-km  segment 
are  shown  in  Figure  1.9.  The  strength  of  turbulence  has  increased  by  more  than  an  order  of 
magnitude,  which  exposes  a  consistent  power-law  regime  between  K=1  and  5  rad/m.  Figures  1.10 
and  1.11  show  the  least-squares  fits  from  sliding  two-decade  ranges  as  obtained  from  the  2^®  blocks. 
The  same  general  characteristics  are  observed,  although  the  slope  estimates  derived  from  the  first 
data  segment  are  more  erratic  when  measured  over  the  two-decade  range. 

In  a  power-law  environment  it  is  desirable  to  isolate  a  spatial-frequency  range  within  a  uniform 
portion  of  the  power-law  region  for  parametric  study.  This  procedure  is  compromised  by  the  lack 
of  clean  power-law  definition  in  the  initial  50  km  of  the  data  set.  Thus,  in  Figure  1.12  the  same 
block  was  analyzed  with  the  two-decade  window  split  into  upper  (red)  and  lower  (blue)  one-decade 
segments.  The  selected  range  for  power-law  extraction  is  K=0.5  to  5  rad/m.  The  two-decade 
fit  (green)  is  replotted  for  reference.  Prom  a  physics  viewpoint,  one  might  expect  a  smooth  low- 
frequency  reservoir  to  feed  the  developing  turbulence.  The  average  spectral  characteristics  shown 
in  Figure  1.13  suggest  this  behavior.  Thus,  the  upper  curve  captures  some  of  the  steepening  while 
the  lower  curve  and  the  overall  average  capture  the  shallower  central  segment  of  the  spectrum.  The 
objective  of  obtaining  a  single  power-law  fit  is  best  met  by  the  shallower  power  law. 

For  completeness.  Figure  1.14  shows  the  average  PSD  for  one  segment  within  the  second  50- 
km  portion  of  the  data.  Here  the  power-law  parameters  are  consistent  and  the  index  is  near  the 
Kolmogorov  value.  We  note,  however,  that  the  lowest  frequencies  depart  from  the  Kolmogorov 
limit  by  first  developing  a  shallow  ledge  and  then  steepening.  The  contamination  of  the  data  at 
low  spatial  wave  numbers  compromises  any  phenomenological  interpretation  of  this  result. 
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ABLE  ACE  95135-04.  nbiock  in  range  4096  to  65536 


Figure  1.8:  Power-law  fit  parameters  (top  slope/bottom  intercept)  derived  from  one-decade  sliding 
window  using  average  spectra  with  different  block  sizes  applied  to  first  50-km  segment. 


K  (rad/m),  (log)center  of  sliding  decade  window 


Figure  1.9:  Power-law  fit  parameters  derived  from  one-decade  sliding  window  using  average  spectra 
with  different  block  sizes  applied  to  second  50-km  segment. 
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ABLE  ACE  95135-04,  nblock  =  32768,  psd  averaged  over  8  blocks  (DispMultWin.m) 


Figure  1.10:  Power-law  fit  parameters  derived  from  two-decade  sliding  window  using  spectra  from 
eight  2^^  sample  blocks  within  the  first  50-km  segment. 


ABLE  ACE  95135-04,  nblock  =  32768,  psd  averaged  over  8  blocks  (DispMultWin.m) 


Figure  1.11:  Power-law  fit  parameters  derived  from  two-decade  sliding  window  using  spectra  from 
eight  2^^  sampleblocks  within  the  second  50-km  segment. 
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Compare  different  tit  ranges  -  smoothed(FitterProcess.m.  fi  le  95135-04 


Figure  1.12:  Power-law  fit  parameters  derived  from  different  spectral  windows  using  spectra  from 
individual  2^^blocks  over  the  entire  data  set. 


ABLE  ACE  951 35-04,  nblock  =  65536,  psd  avg  over  8  blocks  (CompReport.m) 


Figure  1.13:  PSD  estimate  from  the  first  50-km  segment  derived  by  averaging  a  typical  set  of  eight 
PSD  estimates  from  2^®  blocks. 
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ABLE  ACE  95135-04,  nblock  =  65536,  psd  avg  over  8  blocks  (CompReport.m) 


Figure  1.14:  PSD  estimate  from  the  second  50-km  segment  derived  by  averaging  a  typical  set  of 
eight  PSD  estimates  from  2^®  blocks. 
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1.3.2 


Variability  of  power-law  parameters  from  fixed  wave  number  ranges 


Blocks  2^®.  2^^.  2^®  and  2^®  -  unsmoothed  (FilterProcess.m,  file  95135-04) 


Figure  1.15:  Power-law  fit  and  variance  parameters  derived  from  fixed  two-decade  window  applied 
to  contiguous  blocks  over  the  entire  data  set. 

Figure  1.15  shows  single-block  power- law  fit  parameters  derived  from  contiguous  segments  over 
the  whole  data  set.  The  superimposed  curves  are  for  different  block  sizes  and  thereby  show  in¬ 
creasing  levels  of  detail  in  the  parameter  estimates.  The  average  spectral  characteristics  derived 
from  a  fixed  spatial  frequency  range  are  not  affected  by  block  size  because  the  largest  scale  used  for 
parameter  extraction  is  captured  by  the  smallest  block.  The  lower  frame  shows  the  second  moment 
of  the  fluctuations,  which  is  equal  to  the  integral  of  the  PSD.  The  second-moment  data  show  an 
increase  with  block  size  and  a  higher  overall  variability  because  of  the  mix  of  structure  involved. 

As  noted  in  Section  1.2.3,  the  Haar  wavelet  provides  an  intrisically  smoothed  version  of  the  PSD 
estimator  that  faithfully  follows  the  PSD  estimate  for  a  uniform  power-law  segment.  Figures  1.16 
and  1.17  show  a  comparison  of  the  standard  and  wavelet-scale  estimators  at  two  segment  lengths. 
In  general,  the  wavelet-scale-spectrum-based  parameters  are  smoothed  versions  of  the  periodogram- 
spectrum-based  estimators;  however,  in  the  weak  turbulence  regime  there  is  a  small  disparity  in 
the  structure  constant  estimate.  This  may  be  due  to  logarithmic  sampling  effects  as  discussed  in 
Section  1.2.3.  For  the  db5  wavelet  as  shown  in  Figure  1.18,  however,  there  is  a  systematic  bias 
consistent  with  our  observations  of  simulated  data. 
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ABLE  ACE  95135-04,  fitRange=[0.05  0.50],  nufnblocks=4  (wavepsdvarS.m) 


Read  advances  by  1 .31  km  each  point;  Segment  size  =  1 .31  km 


Figure  1.16:  Comparison  of  wavelet-scale-spectrum-based  estimator  to  the  periodogram-based  es¬ 
timator  for  Haar  wavelet  at  1.1-km  resolution. 


ABLE  ACE  95135-04,  frtRange=[0.05  0.50],  numblocks=:2  (wavepsdvat3.m) 


Read  advances  by  0,66  km  each  point;  Segment  size  =  0.66  km 


Figure  1.17:  Comparison  of  wavelet-scale-spectrum-based  estimator  to  the  periodogram-based  es¬ 
timator  for  Haar  wavelet  at  0.6-km  resolution. 
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ABLE  ACE  95135-04,  fitRange=[0.05  0.50],  numbIocks=4  (wavepsdvarS.m) 


Read  advances  by  1.31  km  each  point;  Segment  size  =  1.31  km 


Figure  1.18:  Comparison  of  wavelet-scale-spectrum-based  estimator  using  db5  wavelet  to  the 
periodogram-based  estimator  from  Figure  1.16. 
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1.4  Analysis  of  Bahrain  Data 

To  demonstrate  the  robustness  of  our  analysis  procedure,  we  have  applied  it  to  two  additional  data 
sets  that  have  been  analyzed  in  detail  by  Phillips  Laboratory.  A  full-resolution  PSD  from  the  first 
data  set  is  shown  in  Figure  1.19.  In  this  data  set  there  is  a  high  spatial  wave  number  feature 
that  we  believe  to  be  spurious.  Thus,  the  spectral  analysis  window  was  set  up  over  a  narrower 
window,  as  indicated  in  Figure  1.19.  Following  the  same  procedure  that  was  used  for  the  95135- 
04  data,  the  power-law  index  and  strength  of  turbulence  parameters  were  estimated  using  both 
PSD  and  Haar- wavelet  scale  spectra.  The  results  for  1.2-km  and  0.6-km  segments  are  shown  in 
Figures  1.20  and  1.21,  respectively.  There  are  differences  in  detail  between  these  summary  data 
and  the  corresponding  summary  analysis  performed  by  Phillips  Laboratory.  We  attribute  these 
differences  to  the  fact  that  overlapping  data  segments  were  used  by  Phillips  Laboratory.  Also, 
Phillips  Laboratory  used  16  block  averages  with  linear  trend  removal.  As  already  noted,  to  properly 
interpret  the  parameter  fluctuations,  we  have  avoided  operations  that  remove  any  structure.  Part 
2  discusses  our  procedures  for  characterizing  the  noise  structure. 


ABLE  ACE  C97206-03,  nblock  =  1 6384,  psd  avg  over  4  blocks  (CompReport.m) 


Figure  1.19;  High-resolution  PSD  from  Bahrain  data  set  97206-3. 

Figures  1.22,  1.23,  and  1.24  show  a  representative  PSD  and  the  parameters  that  result  from 
applying  the  parameter  estimation  procediue  to  data  set  97206-15.  Here  the  turbulence  is  somewhat 
weaker  and  the  limitations  of  the  restricted  spectral  window  that  was  available  are  influencing  the 
parameter  estimation.  There  are  significant  difierences  between  the  PSD  and  Haar-wavelet-scale- 
spectra  estimates.  The  logarithmic  sampling  associated  with  the  wavelet  scale  spectrum  produces 
a  sparse  sampling  of  the  available  spectral  window;  moreover,  there  are  significant  departures  from 
power-law  behavior.  This  is  a  topic  for  further  investigation,  but  for  now  we  accept  the  biased 
wavelet  estimate  as  the  reference  for  analysis  in  Part  2. 
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ABLE  ACE  c97206a03,  trtRange-[0.05  0.50],  numblocks=4  (wav6psdvar3.m) 


Read  advances  by  1 .16  km  each  point;  Segment  size  =  1 .16  km 


Figure  1.20:  PSD  and  wavelet-scale  parameter  estimates  from  97206-3  data  set  at  1.2-km  segmen¬ 
tation  length. 


ABLE  ACE  c97206a03,  fitRange=[0.05  0.50],  numblocks=2  (wavepsdvarS.m) 


Read  advances  by  0.58  km  each  point;  Segment  size  =  0.58  km 


Figure  1.21:  PSD  and  wavelet-scale  parameter  estimates  from  97206-3  data  set  at  0.6-km  segmen¬ 
tation  length. 
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ABLE  ACE  C97206-1 5,  nblock  =  16384,  psd  avg  over  4  blocks  (CompReport.m) 
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Figure  1.22:  High  resolution  PSD  from  Bahrain  data  set  97206-15. 


ABLE  ACE  c97206a15,  fitRange=[0.05  0.50],  nunnblocte=4  (wavepsdvarS.m) 


Read  advances  by  1.16  km  each  point;  Segment  size  s:  1.16  km 


Figure  1.23:  PSD  and  wavelet-scale  parameter  estimates  from  97206-15  data  set  at  1.2-km  segmen¬ 
tation  length. 


ABLE  ACE  c97206a15.  fitRange=[0.05  0.50],  numblocks=2  (wavepsdvarS.m) 


Read  advances  by  0.58  km  each  point;  Segment  size  =  0.58  km 


Figure  1.24:  PSD  and  wavelet-scale  parameter  estimates  from  97206-15  data  set  at  0.6-km  segmen¬ 
tation  length. 
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1.5  Summary 

To  summarize  the  results  of  the  preliminary  data  analysis  and  parameter  estimation,  the  parameter 
estimates  from  each  file  at  the  highest  resolution  of  approximately  0.6  m  are  presented  in  histogram 
form.  Figure  1.25  corresponds  to  to  the  data  presented  in  graphical  form  in  Figure  1.17.  Several 
features  should  be  noted.  First,  the  wavelet-scale-spectra  parameters  show  less  fluctuation.  Sec¬ 
ondly,  the  strong  and  the  weak  turbulence  produce  somewhat  distinct  populations  in  the  turbulent 
strength  histograms — one  peaking  below  -18  log  units  and  the  other  peaking  between  -17  and  -16 
log  units.  By  comparison,  the  wavelet  slope  distribution  is  more  uniform.  This  is  an  encouraging 
result  in  that  one  expects  the  slope  parameter  to  be  intrinsic  to  the  turbulence  process.  Departures 
from  the  Kolmogorov  limit  would  occur  only  because  the  data  window  did  not  encompass  mea¬ 
surable  turbulence  in  its  inertial  regime.  Strength,  by  comparison,  should  reflect  how  strongly  the 
turbulence  is  being  driven.  These  observations  apply  to  the  expectation  values  of  the  parameters. 
The  fluctuation  process  will  be  investigated  in  Part  2. 

Similarly,  the  data  presented  in  Figures  1.26  and  1.27  correspond,  respectively,  to  the  graphical 
presentations  in  Figures  1.21  and  1.24.  These  results  show  the  same  smoothing  effect  of  the  wavelet- 
scale-spectrum  estimator,  but  the  distributions  are  more  erratic,  reflecting  the  limited  spectral  range 
discussed  above  and  the  generally  weaker  turbulence  levels  encountered. 
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Figure  1.25:  Histogram  presentation  of  data  from  Figure  1.17. 
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Histogram  lor  segment  =  4  x  Ng=:  16  k 
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Figure  1.26:  Histogram  presentation  of  data  from  Figure  1.21. 
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Figure  1.27:  Histogram  presentation  of  data  from  Figure  1.24. 
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2  Scale  spectral  analysis 


The  basic  scaling  theory  of  atmospheric  turbulence  tells  us  that,  locally,  the  power  spectral  density 
of  the  temperature  should  have  power  law  behavior  over  some  range  of  scales.  The  aim  of  this  study 
is  to  estimate  this  local  power  law  behavior  from  the  aerothermal  data,  allowing  for  variations  over 
longer  scales  that  are  on  the  order  of  kilometers. 

In  Section  2.1  we  introduce  the  data  and  note  that  they  are  corrupted  by  the  measurement 
process.  In  Section  2.2  we  introduce  the  scale  spectral  decomposition  of  the  data,  or  of  segments  of 
it,  and  note  qualitatively  its  power  law  form.  We  discuss  different  segmentations  of  the  data  and  the 
way  a  power  law  is  fitted  to  the  scale  spectra.  In  Section  2.3  we  discuss  how  to  assess  and  filter  out 
the  effects  of  segmentation  and  sampling  variability  in  the  scale  spectra.  We  begin  in  Section  2.3.1 
with  a  simulation  of  the  temperature  data  as  a  Gaussian  local  power  law  process  with  the  slope 
and  log  intercept  obtained  from  the  actual  data  and  then  smoothed  (Figure  2.7).  The  estimation 
of  the  local  power  for  the  simulated  data  gives  us  a  quantitative  measure  of  fluctuations  due  to 
sampling  and  segmentation.  In  Section  2.3.2  we  carry  out  a  variogram  analysis  of  the  estimated 
slope  and  log  intercept  of  the  actual  data.  We  assume  that  the  variograms  have  a  simple  structure 
that  depends  on  three  parameters,  and  we  estimate  them.  In  Section  2.3.3  we  use  the  estimated 
variogram  model  to  filter  out  the  sampling  error  in  the  slope  and  log  intercept  estimates.  We  find 
(Figures  2.14,  2.15  and  2.16)  that  this  filtering  process  removes,  to  a  large  extent,  segmentation 
effects  over  a  large  range  of  segmentation  sizes.  This  filtering  is  a  simple  but  important  step  in  our 
analysis  of  the  aerothermal  data. 

In  Section  2.4  we  show  how  to  simulate  faithfully  the  temperature  data  as  a  local  power  law 
process.  In  Section  2.4.1  we  generate  and  analyze  simulations  which  are  locally  Gaussian  power-law 
processes  over  the  full  range  of  scales  available  in  the  actual  data.  The  local  power  law  that  we 
use  in  the  simulations  is  estimated  one  from  the  actual  data,  after  filtering  (as  in  Figure  2.14).  In 
Section  2.4.2  we  use  replication  to  assess  the  standard  deviation  of  the  sampling  error.  In  Section 
2.4.3  we  examine  the  effect  on  the  estimation  of  the  slope  and  log  intercept  of  the  range  of  scales  of 
the  spectrum  that  is  used.  This  indicates  that  over  a  suitable  range  of  scales  a  local  power  law  model 
fits  the  data  well.  In  Section  2.4.4  we  simulate  the  temperature  data  with  a  Gaussian  local  power 
law  process  in  which  the  slope  and  log  intercept  are  themselves  stochastic  processes  with  statistical 
properties  determined  by  the  actual  data,  as  in  Section  2.3.2  and  the  smooth  background  of  Figure 
2.7.  In  Section  2.4.5  we  use  replication  to  estimate  the  standard  deviation  of  the  sampling  error  in 
the  slope  and  log  intercept  process.  We  find  that  it  is  nearly  twice  as  large  as  the  corresponding 
one  of  Section  2.4.2  where  the  simulations  use  the  actual  estimates  of  the  slope  and  log  intercept. 

In  Section  2.5  we  summarize  our  analysis  and  conclusions.  Our  main  result  is  that  the  atmo¬ 
spheric  temperature  data  that  we  studied  can  be  modeled  very  well  by  a  local  power  law  process. 
The  slope  and  log  intercept  of  the  local  power  law  spectrum  can  be  estimated  accurately  and  are 
essentially  independent  of  the  segmentation  used. 

2.1  The  data 

In  Figure  2.1  (top)  we  show  the  temperature  data,  which  is  composed  of  5,636,096  data  points 
and  is  plotted  subsampled  1/50.  In  Figure  2.1  (bottom)  we  show  the  consecutive  differences  of  the 
temperature,  subsampled  1/50.  It  is  clear  from  the  difference  figure  that  there  is  a  very  pronounced 
periodic  component  in  the  data,  with  period  of  about  4  km  (20  sec).  In  Figure  2.2  we  show  two 
2km  segments  of  the  data,  one  (top)  beginning  at  the  2-km  point  and  the  other  (bottom)  beginning 
at  the  108-km  point.  Figure  2.3  is  the  same  as  Figure  2.2  except  that  temperature  differences  are 
shown.  Note  the  smaller-amplitude,  regular  periodic  pattern  with  period  of  about  60  m  and  the 
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Data;  subsampled  1/50 


Differences;  subsampled  1/50 


Figure  2.1:  Top  frame:  the  raw  temperature  data  over  the  first  84  km  out  of  a  112  km  record,  sub¬ 
sampled  1/50.  The  spatial  resolution  of  the  data  is  approximately  2  cm.  Bottom  frame:  successive 
temperature  differences  for  the  same  84km  data  set,  subsampled  1  /50.  Note  the  systematic  burst 
in  the  tumperature  differences,  with  a  period  of  about  4  km 


larger-amplitude  pattern  of  about  500  m,  which  is  the  blown-up  view  of  the  4-km  blip  in  Figure 
2.1.  This  larger-amplitude  pattern  has  an  internal  periodic  structure  of  about  60  m. 

It  is  clear  from  Figures  2.1,  2.2  and  2.3  that  the  temperature  data  are  not  a  stationary  time 
series  and  do  not  have  stationary  increments  (differences).  However,  the  increments  have  a  more 
regular  structure  and  the  embedded  periodic  components  are  more  visible  in  the  difference  plot. 
It  is  also  clear  that  the  visible  periodic  structures  in  the  data  cannot  be  attributed  to  turbulence, 
or  larger  coherent  structures  in  the  atmosphere,  because  of  the  precise  periodicity  over  the  full 
112-km  length  of  the  data  set.  The  estimation  of  the  effects  of  atmospheric  turbulence  must  seek 
to  suppress  these  spurious  features. 
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Signal;  subsampled  1/50 


Figure  2.2:  Top  frame:  the  temperature  data  between  2  km  and  4  km,  subsampled  1/50.  The 
segment  between  2.4  km  and  2.8  km  corresponds  to  a  burst.  Bottom  frame:  the  temperature  data 
between  108km  and  110km.  There  is  a  burst  between  109  km  and  109,4  km  but  it  is  not  so  visible 
in  this  plot. 
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Figure  2.3:  The  top  and  bottom  frames  here  axe  successive  temperature  differences  for  the  two 
segments  in  Figure  2.2.  Note  that  the  internal  structure  of  the  burst  has  a  repeating  pattern  with 
a  50-  to  60-m  period.  There  is,  in  addition,  a  lower  amplitude,  50-  to  60-m  periodic  pattern  in  the 
rest  of  the  data. 
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2.2  The  wavelet  analysis 

Let  X  =  (ao(l),ao(2),...,ao(2^))  denote  the  temperature  data.  It  is  convenient  to  work  with  data 
vectors  whose  length  is  a  power  of  two  because  computations  are  faster  for  them.  In  our  case  M  is 
between  22  and  23  and  we  shall  work  typically  with  vectors  of  length  2^^  =  4194304,  corresponing 
to  the  first  84-km  section  of  the  112-km  data  set.  We  want  to  carry  out  a  spectral  analysis  of  X 
and  to  fit  the  estimated  power  spectral  density  to  a  power  law.  For  this  purpose  scale  spectra, 
rather  than  Fourier  spectra,  are  a  more  fiexible  tool  and  we  will  use  them  here  [1]. 

From  the  signal,  we  construct  successively  its  wavelet  coefficients  with  respect  to  the  Haar  basis 
as  follows.  Let 


ai(n)  =  •^(oo(2n)  +  ao(2n  -  1)) 

di(n)  =  -^(ao(2n)  —  ao(2n  —  1))  ,  for  n  =  1,2..., 2^“^ 


(2.1) 

(2.2) 


be  the  smoothed  signal  and  its  fiuctuation,  or  detail,  at  the  finest  scale.  Note  that  the  detail  vector 
di  contains  every  other  successive  difference  of  the  data.  This  process  of  averaging  and  differencing 
can  be  continued  by  defining 


and  in  general 


02  (n) 

=  +  ai{2n  -  1)) 

(2.3) 

d2{n) 

=  ■^(ai(2n)  —  ai(2n  —  1))  ,  for  n  =  1,2..., 2^“^ 

(2.4) 

ak{n)  = 

^(afc_i(2n)  +  afc_i(2n  -  1)) 

(2.5) 

4(n)  = 

-^(afc_i(2n)  -  ak-i{2n  -  1))  ,  forn  =  1, 2...,  2^~^ 

(2.6) 

for  k  =  1, ...,  M.  The  data  vector  X  can  then  be  reconstructed  from  qm,  dxj,  du-i,  •••,  di  since  from 
(2.1)  we  have 


“0  (2n)  =  ^  (oi  (n)  +  di  (n) ) 

ao(2n  —  1)  =  -^(ai(n)  —  di(n))  ,  for  n  =  1,2...,  2^“^ 

v2 


(2.7) 

(2.8) 


and  now  ai  can  be  replaced  by  sums  and  differences  of  a2  and  ^2?  etc. 

In  Figure  2.4  we  show  the  first  21  km  of  the  temperature  data  (M  =  20)  along  with 
dyjdg,  dll,  ®12*  Note  the  similaxity  between  di,  and  to  some  extent  d2,  and  the  successive  differences 
in  Figure  2.1,  bottom.  The  detail  coefficients  carry  information  about  the  data  on  longer  scales, 
as  k  increases.  For  example,  dj  shows  successive  differences  of  the  temperature  over  a  distance  of 
2.56  m,  after  averaging  over  successive  segments  of  length  1,28  m. 

The  scale  spectrum  of  X,  relative  to  the  Haar  wavelet  basis,  is  the  sequence  Sj  defined  by 


2M-j 


Sj  =  -2MZj  >  j  =  1,2,...,M 


(2.9) 


n=l 
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Figure  2.4:  Haar  wavelet  coeflScients  for  the  first  21  km  of  the  data.  The  top  left  firame  is  the 
temperature  data.  The  one  below  it  is  the  di  detail  coefficients  (successive  differences  as  in  Fig. 
2.1).  The  third  from  the  top  is  the  and  the  bottom  the  d^  detail  cofficients.  Note  that  the  4-km 
burst  is  not  visible  any  more  at  the  ds  level,  after  four  successive  averagings  of  the  data.  The  top 
frame  on  the  right  is  the  ai2  coefficients  and  below  it  dy,  dg,  dn. 

It  is  easily  verified  from  the  definitions  above  that  the  norm  of  the  data  vector  X  can  be  written 
as 

2^  M 

^(ao(n))2  =  (am)"  +  E  ’  (2-10) 

n=l  j=l 

which  is  another  way  of  expressing  the  orthogonality  of  the  decomposition  of  X  into  gm  and  the 
dj,  j  1,...,JV/. 

In  Figure  2.5  we  show  log-log  plots  of  scale  spectra  of  nonoverlapping  segments  of  the  data,  of 
length  655  m.  Each  plot  contains  scale  spectra  over  15  scales,  from  4  cm  to  655  m.  Note  that  the 
scale  spectra  show  a  distinct  departure  from  power  law  behavior  for  scales  below  one  meter.  On 
the  longer  scale  side,  power  law  behavior  may  be  considered  in  the  range  of  2.5  m  to  80  m,  which 
corresponds  to  the  detail  coefficients  dy  to  du- 

There  are  several  important  issues  that  must  be  addressed  at  this  point,  before  we  can  continue 
with  the  scale  spectral  analysis. 

•  The  temperature  data  are  not  a  stationary  time  series.  Therefore,  computing  scale  spectra 
over  long  data  segments  that  are  not  stationary,  or  approximately  stationary,  gives  Sj's  that 
are  hard  to  interpret  and  probably  meaningless. 

•  The  noise  bursts  in  the  data  enter  only  in  some  of  the  detail  coefficents,  as  can  be  seen  from 
Figure  2.4.  By  restricting  power  law  fitting  to  scale  spectra  from  dy  to  di2  we  eliminate  most 
of  the  noise. 

•  We  want  to  fit  a  power  law  model  to  the  scale  spectra,  like  the  ones  shown  in  Figure  2.5,  but 
even  though  the  average  slope  of  the  log-log  scale  spectra  over  several  segments  is  close  to 
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spectrum.  Section  lenght  655  m 


Figure  2.5:  Scale  spectra  of  655  m  nonoverlapping  segments  of  the  data  (2^^  points  per  segment) 
obtained  from  the  Haar  wavelet  decomposition.  The  dashed  line  has  slope  —5/3  as  in  Kolmogorov 
spectra.  The  solid  line  is  the  average  over  the  scale  spectra  of  the  different  segments.  After  the 
averaging  is  done  the  scale  spectra  are  plotted  in  log-log  format. 


—5/3,  as  the  scaling  theory  of  turbulence  would  predict,  there  is  a  lot  of  variability,  A  local 
power  law  model  is  hkely  to  fit  the  data  much  better. 

We  must  therefore  deal  with  a  larger  class  of  models,  local  power  law  processes,  and  find  effective 
ways  to  estimate  them.  We  will  base  om:  analysis  on  two  criteria.  One  is  that  the  estimation  of 
parameters  should  not  depend  on  how  the  data  are  segmented.  The  second  is  that  the  variations  of 
the  slope  and  intercept  for  the  log  scale  spectra  must  be  on  length  scales  that  are  much  longer  than 
those  used  to  compute  the  spectra  themselves.  We  will  show  with  our  analysis  that  both  criteria 
can  be  met. 

We  begin  by  selecting  four  segmentations  and  estimating  slopes  and  intercepts  for  the  scale 
spectra  obtained  for  each  segment,  as  shown  in  Figure  2.6.  The  four  segmentations  are: 

160m  (2^^  points),  327  m  (2^^),  655  m,  (2^^),  1.31  km  (2^®).  (2,11) 

There  are  512,  256,  128,  64  nonoverlapping  segments,  respectively,  in  each  case.  Within  each 
segment  we  use  only  scale  spectra  from  dy  to  di2,  corresponding  to  length  scales  between  2.5  m  and 
80  m.  We  limit  the  longer  scales  to  80  m,  well  below  the  shortest  segmentation  length  which  is  160 
m.  We  limit  the  shorter  scales  to  2.5  m  because,  as  can  be  seen  from  Figure  2.5,  below  that  scale 
the  spectrum  changes  and  power  law  fit  is  not  appropriate.  We  return  to  the  issue  of  selection  of 
a  range  where  the  power  law  model  is  appropriate  in  Section  2.4.3.  (Figures  2.28,  2,29,  and  2.30). 
The  fitting  is  done  by  least  squares 


log Sj  «  log C  -  plog(|^)  ,  j  =  7, 8, 12 


(2.12) 


for  each  segment,  and  the  tmits  are  adjusted  so  that  the  horizontal  axis  in  Figure  2.6  is  in  kilometers. 
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Slop©  process.  Segmentation  160  m,  328  m,  655  m  and  1 .3km. 


Intercept  process.  Segmentation  160  m,  328  m,  655  m  and  1 .3km. 


Figure  2,6:  Top  frame:  slopes  of  log-log  scale  spectra  from  wavelet  decompositions  based  on  four 
different,  segment  lengths.  At  the  finest  resolution  the  scale  spectra  are  calculated  over  nonover¬ 
lapping  segments  160-m  long.  This  is  the  dotted  line  that  has  the  largest  variations.  At  the  next 
resolution  the  nonoverlapping  segments  are  327-m  long  (also  plotted  with  a  dotted  line).  The  other 
two  resolutions  are  655m  and  1.31  km  and  are  shown  with  dashed  and  solid  lines,  respectively.  The 
bottom  frame  is  the  same  as  the  top,  but  now  for  the  log  intercept  of  the  scale  spectra. 


We  see  from  Figure  2,6  that  the  estimated  slopes,  p,  and  the  log  intercepts,  log  (7,  vary  consid¬ 
erably  over  the  80-km  data  set.  They  also  depend  on  the  segmentation,  with  the  finer  one  having 
larger  fluctuations  (dotted  lines).  We  must  next  develop  a  method  for  removing  this  dependence. 

We  can  extract  from  the  slope  and  log  intercept  process  of  Figure  2.6  the  smooth,  background 
variations.  This  can  be  done  with  a  simple  moving  average  or  with  a  three-term  matching  pursuit 
[2]  and  the  result  is  shown  in  Figure  2,7.  Note  that  the  smooth  variations  do  not  depend  very  much 
on  the  segmentation. 

2.3  Filtering  of  the  parameter  processes 

In  this  section  we  carry  out  a  filtering  of  the  slope  and  intercept  processes  of  Figure  2.6,  which 
approximately  removes  segmentation  effects,  that  is,  the  dependence  of  the  slope  and  intercept  on 
the  length  of  segment  of  the  data  that  are  used  to  estimate  them.  To  motivate  the  filtering,  we 
first  simulate  numerically  temperature  data  with  variable  slope  and  intercept  from  Figure  2.7.  This 
simulation  replicates  some  important  properties  of  the  actual  data.  After  the  filtering,  in  Section 
2,4,  we  will  be  able  to  generate  more  faithful  simulations  of  the  data. 

2.3.1  Synthetic  data 

Stationary,  Gaussian  processes  with  power  law  spectra  can  be  simulated  efficiently  by  the  method 
described  fully  in  [3]  and  briefly  in  Appendix  A.  They  can  be  simulated  by  frequency  and  time- 
domain  methods  as  well.  Local  power  law  processes  are,  of  course,  not  stationary  since  the  slope 
and  the  log  intercept  vary.  But  they  vary  slowly,  on  the  scale  of  hundreds  of  meters,  compared  to 
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Estimated  mean  slopes 


Figure  2,7:  Top  and  bottom  frames  are  smooth  or  background  variations  for  the  slope  and  intercept 
processes  of  Figure  2.6.  The  four  hnes  correspond  to  the  four  different  segmentations  (resolutions). 
The  two  dotted  hnes  correspond  to  the  finer  segmentations  and  the  dashed  and  sohd  fines  to  the 
coarser  ones,  as  in  Figure  2,6. 

the  scale  of  0.2  cm  at  which  the  process  is  sampled.  In  Figure  2,8  we  show  one  realization  of  a 
local  power  law  process  simulated  by  varying  the  slope  and  log  intercept  according  to  Figure  2.7, 
discretized  on  a  1-km  scale.  The  84rkm  simulation  shown  contains  iP"  points. 

The  simulated  temperature  and  temperature  difference  processes  of  Figure  2.8  do  not  look  very 
much  like  the  etctual  data  of  Figure  2,1.  The  simulations  of  Section  2.4  (Figure  2.18  and  especially 
Figure  2,19,  as  well  as  Figure  2.31)  axe,  however,  much  more  faithful  replications  of  the  actual  data, 
over  the  selected  range  of  scales,  as  we  explain  below. 

Even  though  the  simulated  temperature  process  does  not  look  like  the  raw  or  actual  data,  the 
slope  and  log  intercept  processes  that  are  obtained  from  it,  for  the  four  different  segmentations 
(2,11)  we  used  for  the  actual  data,  are  quite  reasonable,  as  can  be  seen  from  Figure  2.9.  The 
fiuctuations  of  the  estimated  slope  and  log  intercept  are  smaller  than  the  ones  for  the  actual  data 
in  Figure  2.6.  This  is  to  be  expected  since  we  used  only  the  smooth  background  of  the  actual  slope 
and  log  intercept.  Figure  2.7,  in  the  simulations.  What  is  important,  however,  is  dependence  on  the 
segmentation.  The  sensitivity  of  the  estimated  slope  and  intercept  processes  to  the  segmentation 
corresponds  roughly  to  that  seen  for  the  actual  data. 

We  can  now  state  more  precisely  our  basic  strategy  for  removing  dependence  on  segmentation 
and  identifying  accurately  the  structure  of  the  temperature  process.  We  will  assume  that  it  is 
a  two-scale  local  power  law  process.  This  means  that  it  is  a  Gaussian  process  with  power  low 
spectrum  for  which  the  slope  and  log  intercept  themselves  are  stochastic  processes  that  vary  slowly 
relative  to  the  underlying  one.  We  will  assume,  moreover,  that  the  slow  modulation  processes  that 
characterize  the  slope  and  log  intercept  are  statistically  independent  of  the  power  law  process. 

We  will  use  the  variogram,  or  structure  function,  to  estimate  statistics  of  the  slope  and  log 
intercept  processes.  For  a  time  series  Y  =  (Fi)  of  size  N  the  empirical  variogram  with  lag  j  is 
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Figure  2.8:  Simulation  of  an  84  km  segment  of  the  temperature  process  as  a  Gaussian,  local  power 
law  process  with  slowly  varying  slope  and  intercept  as  obtained  from  the  actual  data  and  shown  in 
Figure  2.7.  The  simulation  is  done  with  the  algorithm  described  in  Appendix  A.  The  top  frame  is 
the  simulated  temperature  and  the  bottom  its  successive  differences. 

defined  by 

with  dependance  on  the  length  of  the  data  vector  N  not  shown.  In  Figure  2.10  we  show  the 
variograms  for  the  fluctuations  of  the  slope  and  log  intercept  process  of  the  simulated  data.  That 
is,  the  difference  between  the  curves  in  Figure  2.9  and  the  backgroimd  in  Figure  2.7.  On  the  left  in 
Figure  2.10  we  use  the  processes  obtained  with  the  finest  segmentation  (160  m)  and  on  the  right 
with  the  next  coarser  (327  m).  It  is  clear  that  the  variograms  are  essentially  flat,  which  indicates 
that  the  estimated  slope  and  intercept  processes  have  trivial  correlation  structure  (white  noise). 
This  is  not  the  case  with  the  actual  data,  as  we  will  see  in  the  next  section. 

In  Figure  2.11  we  show  Fomier  spectra  for  the  fluctuations  of  the  slope  and  intercept  processes 
of  the  simulated  data.  The  spectra  are  very  fluctuating  but  have  no  particular  structure,  supporting 
thus  the  white  noise  behavior  that  was  suggested  by  the  variograms  in  Figure  2.10. 

2.3.2  Measured  data 

For  the  actual  data  the  estimated  slope  and  log  intercept  processes  are  shown  in  Figure  2.6,  for 
each  of  the  four  segmentations.  The  variograms  for  these  processes  are  shown  in  Figure  2.12  for 
the  finest  segmentation,  and  in  Figure  2.13  for  the  coarsest.  They  clearly  have  more  structure 
than  their  analog  for  the  simulated  data  in  Figure  2.10.  There  are  three  important  parameters 
that  emerge,  and  are  estimated,  from  the  variograms  of  Figme  2.12  that  we  shall  use:  the  vertical 
intercept,  the  level  of  the  horizontal  asymptote,  and  the  decay  rate  of  the  exponential  curve  that  is 
fitted  to  the  points  obtained  from  the  data.  For  the  finest  segmentation  (160  m  each,  512  of  them) 
the  estimated  parameters  are  as  follows: 
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Slope  process.  Segmentation  160  m,  328  m.  655  m  and  1 .3km. 


Intercept  process.  Segmentation  160  m,  328  m,  655  m  and  1 .3km. 


Figure  2.9:  Slope  (top  frame)  and  log  intercept  (bottom  frame)  processes  obtained  from  the  sim¬ 
ulated  temperature  data  shown  in  Figure  2.8.  The  four  curves,  two  dotted,  one  dashed  and  one 
solid,  correspond  to  slopes  and  log  intercepts  obtained  from  scale  spectra  with  segments  of  four 
different  lengths,  as  was  done  for  the  actual  data  and  shown  in  Figure  2.6.  The  smooth  background 
of  these  processes  coincides  with  the  one  we  used  in  the  simulations  (Figure  2.7).  The  fluctuations 
about  the  backgroimd  are  not  as  large  as  the  ones  for  the  actual  data,  shown  in  Figmre  2.6. 

•  For  the  slope:  Vertical  intercept  =0.09,  horizontal  asymptote  level  -|-  crj=0.21,  and 
correlation  length  640  m  (four  segment  lengths). 

•  For  the  log  intercept:  Vertical  intercept  =0.13,  horizontal  asymptote  level  a^+a^  =1.83, 
and  correlation  length  480  m  (three  segment  lengths). 

The  vertical  intercept  can  be  used  to  identify  a  white  noise  component  cr^  in  the  estimated 
slope  and  log  intercept  processes.  The  exponential  and  the  asymptote  level  can  be  used  to  identify 
the  intrinsic  variation  of  the  parameters;  we  model  this  process  component  as  being  stationary 
after  the  backgroimd  is  subtracted  out.  We  also  assume  that  the  white  noise  and  the  exponentially 
correlated  fluctuation  processes  axe  independent.  This  is  a  reasonable  hypothesis  because  of  the 
separation  of  scales  in  the  two-scale  power  law  model  that  we  take  as  modeling  the  temperature 
data.  In  the  next  section  we  use  the  information  we  have  gotten  from  the  variogram  to  introduce 
a  model  for  the  slope  and  intercept  process  of  Figure  2.6  and  then  to  Alter  them  so  as  to  remove 
the  effects  of  segmentation. 

2.3.3  Filtering  of  the  slope  and  log  intercept  processes 

We  now  model  the  slope  and  log  intercept  processes  by  a  white  noise  process  plus  one  with  known 
background  and  covariance  and  then  Alter  out  the  white  noise.  Filtering  out  the  white  noise 
component  of  the  estimated  processes  obtained  from  the  actual  data  will  give  us  the  intrinsic  slope 
and  intercept  processes  that  should  not  depend  on  the  segmentation.  This  is  what  we  do  in  this 
section.  The  slope  and  intercept  processes  for  the  actual  data  are  shown  before  the  filtering  in 
Figure  2.6  and  after  the  filtering  in  Figure  2.14.  Note  that  after  the  filtering  aU  four  segmentations 
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Variogram  for  slope;  fine  segmentation.  Variogram  for  slope;  coarse  segmentation. 


km  km 


Figure  2.10:  Vaxiograms  for  the  fluctuations  of  the  estimated  slope  and  log  intercept  processes  of 
the  scale  spectra  of  the  simulated  temperature  data  at  the  two  flnest  segmentations  (left  frames 
are  for  the  the  flnest  segmentation).  Note  that  the  variograms  are  flat,  which  indicates  that  the 
fluctuations  of  the  slope  and  log  intercept  processes  are  essentially  white. 

give  esentially  the  same  result,  the  solid  curve  in  Figure  2,14.  This  is  a  very  strong  indication  that 
the  model  we  have  chosen  and  the  analysis  that  we  have  followed  fit  the  aerothermal  data  very 
well.  Recall,  however,  that  the  correlation  length  of  the  slope  and  log  intercept  processes  were 
approximately  500  m.  Therefore,  when  the  parameters  are  estimated  based  on  the  two  coarsest 
segmentations  (655  m  and  1.3  km),  the  parameter  processes  are  smoothed  somewhat  and  the 
estimates  do  not  faithfully  represent  their  variability.  This  is  verified  by  Figure  2.14,  since  the  two 
dotted  curves,  corresponding  to  the  coarser  resolutions,  are  somewhat  smoothed  versions  of  the 
solid  line  that  corresponds  to  the  two  finest  resolutions.  In  Figures  2.15  and  2.16  we  show  explicitly 
the  action  of  the  filters.  Figure  2,15  corresponds  to  the  slope  parameter  and  Figure  2,16  to  the  log 
intercept  process.  The  four  frames  in  the  figures  correspond  to  the  diff'erent  segmentations  with  the 
finest  at  the  top.  The  dotted  lines  are  the  estimated  parameters  before  filtering  and  the  solid  lines 
after  filtering.  Note  that  for  the  two  finest  resolutions  the  variablity  of  the  parameters  are  captured, 
but  the  sampling  noise  is  large.  The  sampling  noise  is,  however,  removed  by  the  filtering.  For  the 
two  coarsest  resolutions  the  parameter  processes  are  overly  smoothed.  Figure  2.17  is  a  scatterplot 
of  the  filtered  slope  and  log  intercept  processes.  It  shows  that  there  is  positive  correlation  between 
them,  a  fact  that  has  some  physical  basis. 

The  filtering  is  done  as  follows.  Let  Y  =  (T^),  i  =  1,...,1V'  be  the  observed  slope  or  intercept 
process  obtained  with  a  fixed  segmentation.  We  model  it  by  a  process  of  the  form 

Yi  =  Si  +  7ii  =  Si  +  Sl  +  ni,  i  =  1, ...,  N  (2.14) 

where  S  =  {Si)  is  the  background,  which  we  model  as  deterministic.  S'  =  (5^)  is  the  fluctuation 
about  the  background,  taken  to  be  a  stationary  stochastic  process  with  exponential  correlation, 
and  n  =  (n^)  is  a  white  noise  process  with  fixed  noise  level.  The  three  parameters  needed  in 
the  modeling  are  provided  by  the  variogram  analysis  of  the  previous  section.  We  estimate  the 
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Spectrum  for  slope  residual;  segmentation  160  km. 


Figure  2.11:  Fourier  spectra  of  the  fluctuations  of  the  slope  of  the  scale  spectra  of  the  simulated 
temperature  data.  The  top  frame  is  for  the  flnest  segmentation  and  is  computed  by  averaging  the 
Fourier  spectra  of  the  two  halves  of  the  slope  process,  each  containing  256  points.  The  bottom  frame 
is  for  the  next-to-the-finest  segmentation  and  is  a  single  Fourier  spectrum  for  the  slope  process  with 
256  points.  The  spectra  are  very  fluctuating  but  do  not  have  distinguishing  features. 


background  S  by  the  smooth  curves  shown  in  Figure  2.7. 

We  now  construct  a  filter  (matrix)  A  that  transforms  the  data  Y  to  AY  in  such  a  way  that 


f;{||AY-S||2} 


(2.15) 


is  minimized  over  all  matrices  A  that  also  preserve  the  mean  of  S  of  S,  that  is  AS  =  S.  If  cr^  is  the 
white  noise  level  and  Cg  is  the  covariance  matrix  of  S',  then 

A  =  (Cs  +  anI)~^[Cs  +  (8)  S]  (2.16) 


where  the  vector  v  =  (uj)  is  given  by 


S^(C,  +  a2/)-iS 


(2.17) 


Here  Cs,i  is  the  column  of  the  matrix  Cg  and  the  superscript  T  stands  for  transpose.  Filtering 
of  this  kind  is  discussed  in  [4]. 


2.4  Simulated  data  based  on  the  filtered  slope  and  log  intercept  processes 
2.4.1  Replication  with  filtered  slope  and  log  intercept 

The  filtered  slope  and  log  intercept  processes  shown  in  Figure  2.14  (solid  line)  essentially  do  not 
depend  on  the  segmentation  and  are,  therefore,  intrinsic  to  the  temperature  data.  We  can  now 
simulate  this  data  as  a  Gaussian,  local  power  law  process  with  slope  and  log  intercept  varying 
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Variogram  for  slope 


Figure  2.12:  Variograms  for  the  fluctuations,  relative  to  their  smooth  background  shown  in  Figure 
2.7,  of  the  slope  (top)  and  log  intercept  (bottom)  processes  of  the  scale  spectra  of  the  actual  data, 
obtained  with  the  flnest  segmentation  (Figure  2.6).  The  stars  axe  computed  points  and  the  solid 
curve  a  fltted  exponential  with  a  level  asymptote.  The  important  parameters  obtained  from  the 
fitting  are  the  vertical  intercept  (white  noise  component),  the  decay  rate  (correlation  length),  and 
the  level  of  the  asymptote  (variance  of  the  fluctuations) . 

according  to  Figure  2.14.  We  discuss  the  issue  of  how  to  best  simulate  the  temperature  data  in 
Appendix  A. 

A  simulation  of  the  temperature  data  done  this  way  is  shown  in  Figure  2.18.  This  should  be 
contrasted  to  the  simulation  of  Figmre  2.8  where  only  the  background,  Figmre  2.7,  of  the  slope  and 
intercept  process  is  used.  Figures  2.1  and  2.18  have  common  features,  like  a  more  energetic  second 
half  (in  the  region  of  53  to  83  km). 

A  more  striking  similarity  between  the  actual  and  the  simulated  data  emerges  when  they  are 
both  filtered  so  that  only  scales  dj  —  di2  contribute.  This  corresponds  to  the  length  scales  in  the 
range  2.5  to  80  m  that  we  use  in  estimating  slopes  and  log  intercepts  of  the  power  law.  Figure  2.19 
shows  the  results  and  they  are,  indeed,  striking.  In  Figure  2.20  we  illustrate  that  the  simulation 
replicates  the  data  faithfully  only  for  the  targeted  scales.  The  top  plot  shows  the  reconstruction  of 
the  actual  data  from  di  to  d^  only  and  is  very  different  from  the  corresponding  reconstruction  for 
the  simulated  data,  middle  frame.  Note  also  that  the  reconstruction  of  the  actual  data  from  ai2 
only  is  very  different  from  the  corresponding  reconstruction  of  the  simulated  data,  as  can  be  seen 
from  Figure  2.20,  bottom  frame. 

We  can  now  treat  the  simulated  temperatmre  data  of  Figure  2.18  as  we  did  with  the  actual 
data.  In  Figure  2.21  we  show  the  slope  and  log  intercept  processes  obtained  with  the  four  different 
segmentations  of  (2.11).  These  curves  appear  to  be  very  similar  to  the  ones  for  the  actual  data. 
Figure  2.6,  but  there  are  important  differences.  When  we  compute  the  variogram  of  the  fluctuations 
of  the  slope  and  log  intercept  processes,  the  analog  of  Figure  2.12,  we  get  Figure  2.22. 

The  stars  in  Figure  2.22  correspond  to  the  variogram  obtained  when  the  slope  and  log  intercept 
used  to  generate  the  simulation  (Figure  2.14)  are  subtracted  from  the  corresponding  processes 


37 


0.00 


Variogram  for  slop© 


Figure  2,13:  Same  as  Figure  2.12  but  with  the  slope  and  log  intercept  process  obtained  from  the 
coarsest  segmentation. 


estimated  from  the  simulation.  This  gives  the  variogram  of  the  sampling  error.  Note  that  the 
flatness  and  level  of  the  star  curves  in  Figures  2.10  and  2.22  are  close  to  each  other,  for  the  finest 
resolution.  The  crosses  in  Figure  2,22  correspond  to  the  variogram  of  the  flltereted  slope  and  log 
intercept  when  only  the  smooth  background  is  subtracted  out.  The  solid  line  is  as  in  Figure  2.12, 
which  is  the  model  variogram  for  the  sampled  slope  and  log  intercept  processes  from  the  real  data. 
Note  that  the  crosses  are  located  somewhat  below  this  line.  The  reason  for  this  is  that  the  Altering 
of  the  sampled  slope  and  log  intercept,  which  approximately  removes  the  sampling  error,  reduces 
somewhat  the  variability  in  the  slope  and  log  intercept  processes.  This  is  seen  from  Figure  2.23, 
where  we  replot  the  crosses,  and  the  stars  now  correspond  to  the  variogram  of  the  slope  and  log 
intercept  processes  used  in  the  simulation,  which  are  given  by  Figure  2.14. 

We  see  that  only  the  white  noise  component  persists  (the  flat  hues).  The  solid  curves  in  Figure 
2.22  are  the  ones  from  the  actual  data  of  Figure  2.12.  In  Section  2.4,4  we  show  how  to  generate 
simulated  data  that  replicate  the  actual  variograms.  The  reason  we  do  not  get  this  now  is  because 
we  have  treated  the  filtered  slope  and  log  intercept  curves  of  Figure  2.14  as  deterministic  curves, 
rather  than  realizations  of  a  stochastic  process  with  the  structure  of  Section  2,3.3, 

When  we  now  Alter  the  simulated  slope  and  log  intercept  processes  of  Figure  2.21  to  remove 
the  white  noise  part,  as  we  did  in  Section  2.3.3,  we  get  the  curves  in  Figure  2.24.  These  curves  are, 
again,  essentially  independent  of  the  segmentation  and  dijffer  very  little  from  the  ones  in  Figure 
2.14,  which  we  got  from  the  actual  data,  as  can  be  seen  from  Figure  2.25.  We  have,  therefore, 
recovered  faithfully  from  the  simulated  data  the  information  that  we  used  to  generate  it. 

2,4.2  Sampling  error 

The  slope  and  log  intercept  processes  are  obtained  from  the  temperature  data  and  are,  therefore, 
functionals  of  one  realization  of  a  stochastic  process.  How  much  samphng  error  is  there  in  the 
curves  of  Figure  2.14?  Stated  in  another  way,  if  we  generate  several  independent  realizations  of  the 
temperature  data  that  replicate  faithfully  its  local-power-law-based  statistical  properties,  how  much 
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Figure  2.14:  Filtered  slope  (top)  and  log  intercept  (bottom)  processes  of  the  actual  data  for  the 
four  different  segmentations.  Note  that  aU  four  segmentations  of  the  filtered  processes  are  essen¬ 
tially  the  same.  The  sohd  and  dashed  fines  correspond  now  to  the  two  finest  resolutions  and  are 
indistinguishable.  The  dotted  curves  correspond  to  the  two  coarsest  segmentations  and  are  slightly 
smoother  versions  of  the  other  curves.  The  filtering  has  eliminated  the  differences  in  the  slope  and 
log  intercept  processes  that  are  due  to  the  segmentation  (as  in  Figure  2,6).  It  has  also  eliminated 
the  white  noise  component  of  the  slope  and  log  intercept  fluctuations  that  is  due  to  sampling.  The 
effect  of  sampling  alone  is  seen  from  the  simulated  slope  and  log  intercept  process  in  Figure  2.9. 

variation  will  we  see  in  the  slope  and  log  intercept  processes  from  realization  to  realization?  This 
is  the  bootstrap  way  of  estimating  error.  When  we  do  this  with  twenty  independent  simulations 
of  the  temperature  data,  like  the  one  of  Figure  2,18,  we  get  Figure  2.27,  where  the  mean  and  two 
standard  deviations  about  it  are  shown  (making  the  mean  curve  darker).  The  sampling  standard 
deviation  should  be  close  to  the  noise  level  that  we  got  from  the  variogram  analysis  in  Section 
2.3.2  and  that  we  used  in  the  filtering  in  Section  2.3.3.  This  is  in  fact  the  case  for  the  log  intercept 
process  error  but  not  for  the  slope  error,  which  is  smaller  than  what  we  got  from  the  variogram 
analysis.  This  can  be  explained  by  looking  more  carefully  at  the  bias  error  in  the  slope  in  Figure 
2.25, 
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Fitt®ring  of  tho  slop®  procossos 


Figure  2.15:  Filtered  (solid)  and  unfiltered  (dotted)  slopes  for  the  four  segmentations,  with  the 
finest  (160  m)  at  the  top. 


Filtering  of  the  intercept  processes 
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Figure  2.16:  Same  as  Figure  2.15  but  for  intercepts. 
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Scatterplot  of  filtered  slope  and  intercept 


Figure  2.17:  Scatterplot  of  the  filtered  slope  process  versus  the  filtered  log  intercept  process  for  the 
actual  data  (Figure  2.14).  Note  the  positive  correlation  between  the  slope  and  the  log  intercept. 
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Figure  2.18:  Simulated  temperatmre  (top)  and  temperature  difierences  (bottom)  by  a  Gaussian, 
local  power  law  process  with  slope  and  log  intercept  obtained  from  the  actual  data  and  shown  in 
Figure  2.14  (solid  line).  If  only  the  background  (Figure  2.7)  is  used  the  simulation  is  shown  in 
Figure  2.8.  The  differences  between  the  two  simulations  are  more  apparent  in  the  temperature 
differences  (the  bottom  frame  here  and  in  Figure  2.8). 
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Figure  2.19:  The  actual  temperature  data  set  (top  frame)  obtained  by  synthesizing  it  from  the 
detail  coefficients  dj  —  di2  of  its  Haar  decomposition.  It  is  information  from  this  range  of  scales 
that  enters  into  the  slope  and  intercept  processes.  The  bottom  frame  is  the  same  as  the  top  but 
for  the  simulated  data  of  Figure  2.18.  Note  the  striking  similarity  between  the  two  figures. 
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Figure  2.20:  The  temperature  data  reconstructed  with  only  the  d\  —  ds  coefficients  (top  frame). 
The  middle  figure  is  the  same  as  the  top  but  for  the  simulated  data  of  Figure  2.18.  Note  the 
dissimilarity  between  the  two  figures.  The  bottom  frame  is  temperature  data  reconstructed  with 
only  ai2,  for  the  actual  data  (solid)  and  the  simulated  data  (dashed). 
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Slope  process.  Segmentation  160  m,  320  m,  655  m  and  1 .3km. 


Intercept  process.  Segmentation  160  m,  320  m,  655  m  and  1 .3km. 


Figure  2.21:  Slope  (top)  and  log  intercept  (bottom)  processes  for  scale  spectra  of  the  simulated 
process  of  Figure  2.18,  with  four  different  segmentations,  as  with  the  actual  data  and  shown  in 
Figure  2.6.  The  Figure  above  and  Figure  2.6  are  quite  similar. 


Variogram  for  slope  Variogram  for  slope 


km  km 


Figure  2.22:  Varigrams  for  the  fluctuations  of  the  slope  (top)  and  log  intercept  (bottom)  processes 
obtained  from  the  simulated  temperature  data  shown  in  Figure  2.18.  The  stars  are  computed  from 
the  fluctuations  when  the  filtered  background  is  subtracted.  The  solid  curves  are  exactly  the  ones 
fitted  to  the  actual  data  and  shown  in  Figure  2.12.  The  crosses  are  variograms  for  fluctuations 
when  only  the  smooth  backgroimd  (Figure  2.7)  is  subtracted.  The  curves  on  the  left  are  from  data 
at  the  finest,  160  m,  resolution  while  those  on  the  right  are  for  the  next  to  finest,  330  m. 
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Variogram  for  slope 


km 


Figure  2.23:  The  crosses  are  as  in  Figure  2.22  and  the  stars  are  the  empirical  variogram  of  the  actual 
data  that  are  used  in  the  simulation,  that  is,  variograms  for  the  slope  and  intercept  of  Figure  2.14 
after  filtering. 


Predicted  slope  processes  for  the  4  segmentations. 


Figure  2.24:  Filtered  slope  (top)  and  log  intercept  (bottom)  processes  from  the  scale  spectra  of 
the  simulated  temperature  shown  in  Figure  2.18.  The  different  curves  correspond  to  the  different 
segmentations,  as  in  Figure  2.14  for  the  actual  data.  The  similarity  between  the  cmrves  above  and 
the  ones  in  Figure  2.14  show  that  the  simulations  rephcate  faithfully  the  information  from  the 
actual  data  that  is  used  in  generating  the  simulated  data. 
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Error  in  slope  and  intercept  process 


Figure  2.25:  The  top  frame  shows  two  curves,  the  slope  process  for  the  actual  data  (Figure  2.14, 
top)  and  the  slope  process  for  the  simulated  data  (Figure  2.24,  top).  The  second  figure  from  the 
top  is  the  difference  between  these  two  curves.  Note  the  slight  bias  in  the  otherwise  small  error. 
The  bottom  two  frames  are  the  same  as  the  top  two  but  for  the  log  intercept  process. 


Scatterplot  of  filtered  slope  and  intercept 


Figure  2.26:  Scatterplot  for  the  filtered  slope  and  intercept  process  of  Figure  2.24  (simulated  data). 
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Confidence  band  for  slope  estimate 

Or 


Figure  2.27:  Mean  over  twenty  independently  simulated  filtered  slope  (top)  and  log  intercept 
(bottom)  processes.  Each  simulation  of  the  temperature  data  is  like  the  one  of  Figure  2.18.  Each 
filtered  slope  and  log  intercept  process  is  like  Figure  2.24.  Twice  the  pointwise  standard  deviation 
is  shown  (to  scale)  just  above  the  horizontal  axis  of  each  figure.  The  two-standard-deviation  band 
is  also  superimposed  on  the  mean  values  of  the  slope  and  intercept,  making  them  appear  darker. 
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2.4.3  Error  in  parametric  fitting  of  scale  spectrum 

We  explained  in  Section  2.2,  below  (2.11),  that  the  slopes  and  log  intercepts  are  obtained  by  fitting 
the  log  scale  spectra  from  dr  to  du  to  a  line,  in  each  segment.  How  good  is  this  fit  and,  more 
important,  could  we  have  done  better  with  a  wider  or  a  narrower  range  of  scale  spectra?  We 
measure  the  error  in  the  fit  by  the  sum  of  the  squares  of  the  differences  of  the  log  spectra  from  the 
fitted  hne.  For  each  segment  we  then  divide  this  error  by  the  average  over  the  different  segments  of 
the  corresponding  error  for  the  simulated  temperature  data.  This  is  shown  in  the  top  Figure  2.28. 
If  for  the  simulated  data  we  replace  the  numerator  by  the  error  in  the  fit  in  each  segment,  we  then 
get  the  bottom  frame  in  Figure  2.28.  The  fact  that  the  top  and  bottom  frames  in  this  Figure  are 
comparable  shows  that  the  model  is  appropriate  over  the  considered  scales.  That  is,  the  misfit  in 
the  actual  data  (top  frame)  is  comparable  to  the  one  for  simulated  data,  which  is  the  most  we  can 
expect. 

In  Figure  2.29  we  show  what  happens  when  we  use  two  additional  small  scale  spectra  in  the  fit. 
This  takes  us  to  the  right  of  the  change  point  (at  about  60  cm)  in  Figure  2.5  so  we  expect  a  worse 
fit.  This  is  in  fact  what  Figure  2.29  shows,  since  the  relative  error  for  the  actual  data  (top  frame) 
is  substantially  bigger  than  the  simulated  (reference)  relative  error  (bottom  frame). 

The  situation  is  similar  if  we  use  only  the  smooth  background  (Figure  2.7)  and  calculate  the  fit 
in  each  segment  with  values  from  it.  We  then  divide  this  error  with  the  corresponding  one  for  the 
simulated  data  where  only  the  smooth  background  is  used  (Figure  2.8).  We  see  that  for  the  actual 
data  the  relative  error  (top  frame)  is  substantially  bigger  than  the  relative  error  for  the  simulated 
data,  the  reference  error  (bottom  frame).  Hence,  a  local  power  law  model  is  not  a  good  model 
when  we  consider  these  extended  sets  of  scales. 
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Ratio  between  observed  and  mean  synthetic  model  fit 


0  10  20  30  40  50  60  70  80  90 


Figure  2.28:  Top:  The  ratio  of  the  f'  norm  of  the  difference  between  the  log-log  scale  spectrum 
and  the  power-law  fitted  log-log  spectrum  for  the  actual  data  (finest  segmantation,  2^^  points  each) 
over  the  mean  (over  the  512  segments)  of  the  norm  of  the  difference  between  the  scale  spectrum 
and  the  power-law  fitted  spectrum  for  the  simulated  data.  Bottom:  Same  as  for  the  top  frame  with 
only  simulated  data.  The  scale  spectrum  corresponding  to  d-j  to  di2  is  used  in  calculating  the  P 
error  and  in  the  power  law  fit.  This  contains  information  in  the  data  over  length  scales  of  2.5  m  to 
80  m. 


Ratio  between  observed  and  mean  synthetic  model  fit 


[km] 


Figure  2.29:  Same  as  Figure  2.28  except  that  two  additonal  scale  spectral  points  are  used,  corre¬ 
sponding  to  dg  and  dg.  This  contains  information  in  the  0.6  m  to  80  m  range. 
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Ratio  between  observed  and  mean  synthetic  model  fit 


Figure  2,30:  Same  as  in  Figure  2.28  except  that  the  simulations  are  based  on  the  smooth  back¬ 
ground,  as  in  Figure  2.7,  and  the  residual  is  also  computed  relative  to  the  smooth  background.  Note 
that  the  relative  error  for  the  actual  data  (top  frame)  is  considerably  larger  than  the  corresponding 
one  for  the  simulated  data  (bottom  frame).  This  should  be  contrasted  with  Figure  2.28  where  the 
two  error  ratios  are  comparable. 
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2.4.4 


Replication  with  random  slope  and  log  intercept 


Realization  I;  subsampled  1/50 


Figure  2.31:  Two  different  realizations  of  simulated  temperature  data  (top  and  bottom  frames)  in 
which  the  slope  and  log  intercept  are  Gaussian  stochastic  processes  with  mean  as  in  Figure  2.7  and 
fluctuations  determined  by  the  variogram  analysis  of  Section  2.3.2.  The  slope  and  the  log  intercept 
fluctuations  are  generated  from  the  same  white  noise  sequence  that  is  independent  of  the  power 
law  process. 

The  replication  of  the  temperature  data  as  a  local  power  law  process  with  slowly  varying, 
deterministic  slope  and  log  intercept  is  done  in  Section  2.4.1.  The  Altered  slope  and  log  intercept  of 
Figure  2.14  are  used  and  a  typical  reahzation  is  shown  in  Figure  2.18.  Simulations  of  this  sort  are 
much  more  faithful  than  the  ones  based  only  on  the  smooth  background  of  Figure  2.7,  and  shown 
in  Figure  2.8.  However,  even  though  the  simulations  of  Section  2.4.1  rephcate  well  the  input  slope 
and  log  intercept,  as  seen  from  Figure  2.25,  the  output  slope  and  log  intercept  of  Figure  2.24  do 
not  have  the  right  correlation  structure.  This  is  clear  from  the  variogram  analysis  of  the  simulated 
slope  and  log  intercept  fluctuations  shown  in  Figure  2.22.  The  filtering  of  the  processes  removes 
some  of  the  intrinsic  variabihty  of  the  parameter  processes,  which  puts  the  crosses  considerably 
below  the  solid  hues. 

How  can  we  generate  realizations  of  temperature  data  that  produce  slope  and  log  intercept 
processes  with  a  variability  like  the  ones  for  the  actual  data  in  Figure  2.12?  To  this  effect,  we  let 
the  slope  and  log  intercept  be  themselves  stochastic  processes  that  have  the  form  (2.14),  with  the 
background  S  given  as  in  Figure  2.7,  fluctuations  about  the  background  S'  taken  to  be  stationary, 
Gaussian  stochastic  processes  with  exponential  covariance  as  estimated  in  Section  2.3.2,  and  without 
the  white  noise  term  n.  The  fluctuation  processes  for  the  slope  and  log  intercept  are  generated 
from  the  same  white  noise  sequence,  which  we  take  to  be  independent  of  the  underlying  power  law 
process.  This  is  a  crucial  assumption  that  can  be  justified  by  the  wide  separation  scales  between 
variations  of  the  basic  power  law  process  (centimeters  to  tens  of  meters)  and  those  of  its  slope  and 
log  intercept  (hundreds  of  meters  to  kilometers).  We  use  the  same  white  noise  sequence  to  generate 
both  the  slope  and  log  intercept,  since  we  expect  these  processes  to  be  strongly  correlated. 

Simulated  temperature  data  with  random  slope  and  log  intercept  are  shown  in  Figure  2.31. 
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Differences  for  'deterministic'  simulation;  subsampled  1/50 


Figure  2.32:  Absolute  values  of  temperature  differences  for  the  simulation  of  Figure  2.18  with 
deterministic  slope  and  intercept  (top  frame)  and  for  the  top  of  Figure  2,31  (bottom  frame). 


Comparison  with  Figure  2.18,  which  is  simulated  data  that  has  deterministic  slope  and  log  intercept, 
shows  that  the  two  data  sets  are  qualitatively  similar  but  the  fluctuations  in  the  differences  are 
bigger  for  the  random  case,  as  can  be  seen  from  Figure  2.32  where  we  plot  the  absolute  value  of 
temperature  differences  for  the  case  of  deterministic  slope  and  log  intercept  (top  frame)  and  the 
random  case  (bottom  frame).  The  estimated  slope  and  log  intercept  for  the  simulation  of  Figure 
2,31,  for  the  four  segmentations  (2.11),  is  shown  in  Figure  2.33.  The  variograms  of  the  fluctuations 
of  the  slope  and  log  interface  of  Figure  2.33  are  shown  in  Figure  2.34  (stars).  We  see  from  this 
figure  that  the  estimations  are  slightly  biased,  which  is  not  surprising  since  complex  nonlinear 
transformations  are  involved.  The  scatterplot  for  the  estimated  slope  and  log  intercept  is  shown 
in  Figure  2,36  and  it  does  indicate  positive  correlation  as  in  Figure  2.17.  The  filtered  slope  and 
log  intercept  processes  are  shown  in  Figure  2.35.  From  Figure  2,24  we  see  that  the  fluctuations  in 
the  slope  are  comparable  in  the  two  simulations,  while  the  fluctuations  of  the  log  intercepts  are 
stronger  for  the  ‘random’  case  of  Figure  2.31. 

It  is  clear  that  when  the  input  slope  and  log  interface  in  the  simulation  of  a  local  power  law 
process  are  themselves  random,  we  do  not  get  estimated  slopes  and  log  intercepts  that  replicate 
pointwise  the  actual  ones.  So  we  should  not  expect  a  close  match  between  Figures  2.14  and  2.35, 
as  we  have  between  Figure  2.14  and  2.24  (see  Figure  2.37).  We  do,  however,  get  a  close  match  (up 
to  a  small  bias)  of  the  variograms  of  their  fluctuations,  the  stars  and  crosses  in  Figure  2.34.  So  the 
statistics  of  the  slope  and  log  intercept  are  now  well  replicated. 

2,4.5  Intrinsic  variability  of  the  slope  and  log  intercept  process 

Let  us  also  look  at  the  standard  deviation  of  the  parameter  estimates  by  generating  several  (twenty) 
independent  realizations,  as  in  the  previous  section,  and  then  estimating  and  filtering  the  slope  and 
log  intercept  processes.  As  shown  in  Figure  2.38  the  two-standard-deviation  band  is  nearly  twice 
as  big  as  that  obtained  with  simulations  that  have  ‘deterministic’  parameters,  as  in  Figure  2.27, 
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Slope  process.  Segmentation  160  m,  328  m,  655  m  and  1.3km. 


Intercept  process.  Segmentation  160  m,  328  m,  655  m  and  1.3km. 


Figure  2.33:  Slope  (top)  and  log  intercept  (bottom)  processes  estimated  for  the  four  basic  segmen¬ 
tations  (2.11)  iOrom  the  simulated  data  of  Figure  2.31. 

2.5  Summary  and  conclusions 

The  main  premise  driving  our  analysis  of  temperature  data  from  a  turbulent  atmosphere  is  that  it 
is  a  local  power  law  process.  This  means  that  the  power  law  itself-the  power  (slope)  and  the  mul¬ 
tiplicative  constant  (log  intecept)-is  not  a  constant  but  a  slowly  varying  function,  deterministic  or 
random.  We  estimate  the  slope  and  log  intercept  of  the  scale  spectrum  by  appropriately  segmenting 
the  data  and  then  removing  segmentation  effects  by  a  filtering  process.  When  a  local  power  law 
process  is  simulated  with  the  parameters  estimated  from  the  data,  we  get  a  faithful  rephcation  of 
its  basic  statistical  properties.  The  replication  is  sample-wise  faithful  but  statistically  inaccurate 
when  the  slow  variation  is  deterministic  (Section  2.4.1).  The  simulation  is  sample- wise  distinct 
from  the  actual  data  when  the  slope  and  log  intercept  are  taken  as  random  in  the  simulation, 
but  the  statistics  are  faithfully  reproduced  (Section  2.4.4).  This  and  the  removal  of  segmentation 
effects  are  a  very  strong  indication  that  the  model  and  the  analysis  that  we  have  followed  capture 
the  essential  features  of  the  data. 

An  important  aspect  of  the  model  that  we  use  is  separation  of  scales  in  the  variation  of  the 
estimated  parameters  (slope  and  log  intercept)  from  the  imderlying  process  that  generates  the 
power  law  spectra.  This  will  be  the  starting  point  for  a  detailed  theoretical  development  of  the 
methods  that  we  have  introduced  here. 
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Figure  2.34:  Variograms  for  the  slope  and  intercept  fluctuations  from  Figure  2,33,  with  the  finest 
resolution.  The  stars  are  from  the  simulated  data  and  the  crosses  from  the  actual  data  of  Figure 
2.12.  Note  that  the  horizontal  asymptote  for  the  slope  is  under-estimated  while  the  one  for  the  log 
intercept  is  over-estimated.  This  is  due  to  a  small  bias  in  the  nonhnear  estimation  procedure. 


Predicted  slope  processes  for  the  4  segmentations. 
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[km] 

Predicted  intercept  processes  for  the  4  segmentations. 


Figure  2.35:  The  filtered  slope  and  log  intercept  processes  obtained  from  Figure  2.33  as  described 
in  Section  2.3.3,  This  is  comparable  to  Figure  2.14  and  2.24  insofar  as  segmentation  effects  have 
been  removed.  It  is,  however,  another  realization  of  the  slope  and  log  intercept  processes  with  the 
correct  statistical  behavior. 


Scatterplot  of  filtered  slope  and  intercept 


Figure  2.36:  Scatterplot  for  the  estmated  slope  and  log  intercept  processes  of  Figure  2.35  (simulated 
data  with  random  slope  and  intercept).  This  is  like  Figures  2.17  and  2.26.  Note  the  positive 
correlation. 


Difference  between  slope  and  intercept  processes 


Figure  2.37:  This  is  analogous  to  Figure  2.25  but  now  the  slope  and  intercept  for  the  actual  data 
is  compared  to  that  obtained  from  the  simulation  with  random  slope  and  intercept  as  in  Figure 
2.36.  The  top  frame  shows  two  curves,  the  slope  process  for  the  actual  data  (Figure  2.14,  top)  and 
the  slope  process  for  the  simulated  data  (Figure  2.36,  top).  The  second  figure  from  the  top  is  the 
difference  between  these  two  curves.  Note  the  size  of  the  difference  compared  to  that  in  Figure 
2.25.  The  bottom  two  frames  are  the  same  as  the  top  two  but  for  the  log  intercept  process. 
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Figure  2.38:  One  realization  of  the  estimated  slope  (top)  and  intercept  (bottom)  of  the  simulated 
temperature  data  with  random  slope  and  intercept.  The  shaded  region  at  the  bottom  of  each  figure 
represents  two  standard  deviations  about  the  values  used  in  the  simulation,  averaged  over  twenty 
reahzations.  This  Figure  should  be  compared  to  2.27  where  the  same  slope  and  intercept  is  used 
in  all  twenty  realizations.  The  standard  deviation  is  larger  now.  The  two  dashed  curves  in  the  top 
and  bottom  firame  correspond  to  the  slope  and  log  intecept  used  in  the  simulation  while  the  solid 
lines  correspond  to  those  estimated  after  the  simulation. 
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A  Simulation  of  power  law  processes 


In  this  appendix  we  show  how  the  realizations  of  the  local  power  law  model  are  being  generated. 
We  differentiate  between  the  cases  that  the  parameter  sequence  (log-intercept  and  slope  of  the 
power  law)  is  given,  and  the  case  where  we  in  a  first  step  draw  this  from  an  appropriate  stochastic 
model.  With  some  abuse  of  terminology  we  refer  to  the  first  as  deterministic  simulation  and  the 
latter  as  random  simulation.  In  the  next  sections  we  define  the  stochastic  model  for  the  local  power 
law.  Thereafter  we  outline  the  simulation  algorithm. 

The  simulation  algorithm  corresponds  to  the  one  presented  in  [3].  In  [3]  our  main  objective 
was  to  define  a  simulation  algorithm  applicable  for  simulation  a  multiple  dimensions.  Here  we  only 
consider  one-dimensional  simulation,  but  use  parameters  that  vary  over  the  simulation  interval, 
corresponding  to  a  local  power  law.  In  the  random  simulation  these  parameters  are  realizations 
from  appropriate  stochastic  models. 

A.l  Model  of  the  process 

Denote  the  stochastic  process  Y  and  a  particular  realization  for  ?/.  We  sample  the  process  on  a 
regular  grid,  the  grid  values  are  yi  for  1  <  z  <  2^  =  iV.  The  nodes  we  denote  Xi  =  i  Ax. 

Let  a  Gaiissian  intrinsic  power  law  process  Y  be  defined  by 


E[{Yi+k-Yi)^]  =  C\xi+k-Xi\P  (A.1) 

and  by  the  requirement  that  the  increments  are  zero  mean  Gaussian  random  variables.  Hence,  the 
marginal  distribution  of  the  samples  need  not  be  Gaussian.  However,  here  we  shall  condition  our 
realizations  on  yo?  hence  the  samples  will  be  Gaussian  and  have  bounded  variance. 

We  next  generalize  (A.l)  to  a  local  power  law,  the  model  introduced  in  the  paper.  For  a  local 
power  law  process  Y 


E[{Yi+k-Yi?]  =  C{xm)\xi^k-xi\^^^-^  (A.2) 

It  is  from  this  model  we  want  to  draw  realizations.  As  mentioned  above  the  smooth  parameter- 
sequences  C{x)  and  p{x)  are  either  defined  a-priori,  or  they  are  draw  from  a  stochastic  model. 
We  model  the  parameters  as  a  truncated  Gaussian  exponentially  correlated  vector  process.  The 
stochastic  model  for  the  parameters,  p,  is  defined  by 


E[{p{xi+k)  -  p{xi+k))'^ {p{xi)  -  p{xi))]  =  Cexp{-kAx/a)  (A.3) 

P  =  ^{p) 

with  C  being  a  covariance  matrix  and  E  a  truncation  operator.  For  the  truncation  Umits  we  used, 
set  so  as  to  ensure  positivity,  none  values  were  truncated.  Note  that  the  ‘smooth-background’ 
or  mean  of  the  parameters  vary  with  location,  whereas  the  random  residual  is  assmned  to  be 
stationary.  Next  we  consider  simulation  from  model  (A.l).  This  will  serve  to  motivate  our  approach 
for  simulation  from  (A.2). 
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A. 2  Direct  simulation  of  a  stationary  sequence 

From  (A.l)  we  can  easily  construct  the  covariance  matrix  of  the  sequence  {Iq,  •  •  •  ,F}v},  denote  it 
Cy.  Then,  a  sample  from  the  process  can  simply  be  computed  using  the  Cholesky  factor  Gy  of  Cy 
[Cy  =  G.GJ)  by 

y  <—  Gy  •  A/”  +  2/0  (A.4) 

with  Af  being  a  vector  of  independent  unit-variance  Gaussian  samples.  This  will  in  general  be 
computationally  too  costly.  We  seek  therefore  an  alternative  simulation  procedure,  the  procedure 
will  be  hnear  in  the  number  of  nodes.  We  simulate  the  process  by  a  sequential  simulation  procedure. 

A. 3  Sequential  simulation  of  a  stationary  sequence 

The  simulation  in  (A.4)  can  be  implemented  as  a  sequential  simulation  in  the  following  way.  Sim¬ 
ulate  the  samples  2/1  •  •  •  2//V  in  sequence,  each  sample  being  conditioned  on  the  previously  simulated 
samples.  That  is  (using  MATLAB-like  notation) 

2/(0  :  AT)  =  0; 
for  i  =  1  :  N 

y{i)  =  g{i,  1  :  i  -  1)  ■  y{l  :  i  -  1)  +  g{i, i)  N{i); 
endy  =  2/  +  2/o; 

with  g{i, :)  being  the  row  vector  containing  the  appropriate  coeflScients  for  the  conditional  simulation 
and  N  a  vector  of  realizations  of  independent  unit  variance  Gaussian  random  variables.  With  the 
right  ordering  in  the  simulation,  the  correlation  with  previous,  distant  samples  might  be  small.  In 
the  markov  case,  each  sample  need  only  be  conditioned  on  the  previous  sample,  this  one  sort  of 
‘screens’  the  effect  of  the  previous  values.  E  we  let  Gi  denote  the  matrix  that  is  the  identity,  apart 
from  the  i^th  row  that  is  defined  by  g{i^0  :  i),  then 


N 

G,  =  Y[Gi.  (A.5) 

1=1 

We  seek  to  construct  the  simulation  sequence  such  that  each  sample  need  only  be  conditioned  on 
a  few  previous  samples,  or  such  that  g{i^  0  :  i  —  1)  only  has  a  few  dominant  terms. 

Consider  first  a  trimcated  version  of  a  classical  method  for  constructing  Brownian  motion,  see 
for  instance  Section  2,3  [5].  This  will  show  how  to  construct  the  right  ordering  such  that  we  can 
condition  on  a  small  set  of  previous  samples  despite  the  fact  that  we  are  dealing  with  long-range 
processes. 
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%  Classical  method  for  constructing  Brownian  motion 
function  y  =  screen{M,  dx,  yO) 

% 

%  Input 

%  M  :  dyadic  dimension  (#grid  points  is  N  =  2^) 

%  dx  ;  grid-spacing 
%  yO  :  value  of  process  at  origin 
%  Output 

%  y  :  simulated  process;  column  N-vector 

% 

JV  =  2^; 

ranN  =  randn{N  +  1, 1); 
y(l)  =  yO; 
cr  =  ViV  dx; 

y{N  +  1)  =  yO  -h  <7  ranN{N  +  1); 
for  i  =  1  :  M 
j  = 

<7=  VWIM; 

/or  fc=l+j:2j:iV  +  l““j 
y{k)  =  {y{k- j)+y{k-\- j))/2  + a  ranN {k)\ 
end 
end 

y  =  y{2:N  +  l); 

The  simulation  scheme  starts  out  by  initializing  the  boundary  nodes.  The  algorithm  then  proceeds 
by  drawing  the  samples  on  dyadic  subgrids,  creating  a  process  of  successively  finer  resolution. 
The  main  aspect  of  the  procedure  being  that  when  the  grid  is  refined,  the  new  samples  (one  at 
each  midpoint  between  samples  at  the  previous  level)  need  only  be  conditioned  on  the  two  nearest 
neighbors^  both  being  simulated  at  the  previous  level.  Thus,  with  this  ordering  in  the  simulation, 
the  vector  gi  contains  only  two  non-zero  entries. 

Consider  next  simulation  of  firactional  Brownian  motion  p  1.  Then  g(i,0:i-l)  will  in  general 
have  more  than  two  non-zero  elements.  However,  sequential  simulation  with  the  same  sequencing 
of  nodes  and  the  same  conditioning  sets  as  in  the  above  construction  gives  accurate  realizations 
also  when  we  consider  fractional  Brownian  motion.  ‘Accurate’  in  the  sense  of  having  essentially 
the  right  moments,  see  [3].  That  is,  we  approximate  the  ‘Cholesky  factor’  Gi  with  a  factor  having 
only  two  elements  in  its  i^th  row,  in  addition  to  the  diagonal  element.  Additional  elements  does  not 
change  the  statistics  by  much.  The  diagonal  element  and  the  two  off  diagonal  elements  are  chosen 
such  that  the  sample  has  the  right  variance  and  expectation  given  the  two  neighbor  elements. 
From  symmetry  the  off  diagonal  coefficients  equal  1/2.  Hence  the  only  change  in  screen  is  that 
the  expressions  for  the  conditional  standard  deviation  a  is  slightly  changed  and  now  depend  on  the 
values  of  p  and  C,  see  [4]. 

We  next  describe  the  specifics  of  the  algorithm  we  used  to  draw  realizations  from  the  model 
A.2. 
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A.4  Screening  simulation  in  non-stationary  case. 

Assume  first  that  the  paxameter  sequences  are  given.  The  parameter  sequences  characterize  the  log- 
mtercept  and  slope.  These  are  the  estimated  slopes  and  log-intercepts  for  the  finest  resolution, 

egments  of  length  2  ,  for  p-13.  The  total  grid  size  is  2^,  with  M  =  22  in  our  case.  The  algorithm 
consists  of  two  steps: 

Step  I 

De^e  a  ‘super-grid’  of  size  2(^-)  +  1  and  draw  the  samples  on  this  according  to  a  model  as  in 
fi  ^  “Odes  m  the  super-grid  corresponds  to  the  boundary  nodes  for  the  segments  of  length 

and  .St  S  T  m  ^  of  the  slope 

og  m  orcept.  Note  that  the  size  of  the  supergrid  is  so  small  that  we  can  actually  use  the 

ecomposition  (A.4).  We  chose  to  do  the  simulation  on  the  supergrid  sequentially.  The  value  of 

HoSTSTp  on  the  previous  boundary  node  using  the  parameter  values 

(log  intercept  and  slope)  that  corresponds  to  the  section  in  between  them. 

Step  II 

In  each  section  of  the  supergrid  simulate  a  firactional  Brownian  motion  by  the  screening  algorithm. 
The  simulation  is  conditioned  on  the  two  boundary  nodes  (in  the  supergrid). 

Thus  the  function  ret  inning  the  simulated  process  is: 

%  Simulation  of  local  power  law 

function  y  =  SimLpl{M^  dx,  yO,  logint^  slope) 

% 

%  Input 

%  M  :  dyadic  dimension  (#grid  points  is  JV  =  2^) 

%  dx  :  grid-spacing 
%  yO  :  value  of  process  at  origin 

%  logint  :  the  log-intercept  parameter  values;  column  vector 
%  slope  :  the  slope  parameter  values;  column  vector 
%  Output 

%  y  :  simulated  process;  column  N-vector 

% 

AT  =  2^; 
y  =  zeros{N^iy^ 
segments  =  length{slope)] 
k  =  N /  segments] 

P  =  log2{k)] 
fori  =  l:  segments 
zT  =  1  +  ^  (z  -  1); 
z2  =  k  i] 

y{il  :  i2)  =  screen{p,dx,yO,logint{i),slope{i)y, 
yO  =  y{i2); 
end 

iTT  screen  is  defined  exactly  as  above  only  that  the  expressions  for  a  is  shghtly  modified 

using  the  relevant  values  of  logint  and  slope.  moamea 
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Finally,  if  we  want  to  carry  out  a  random  simulation  then  the  vectors  logint  and  slope  are 
sampled  from  the  model  (A.3).  Simulation  from  this  model  is  easily  implemented  since  p  is  a 
Markov  process. 

B  Analysis  of  other  aerothermal  data  sets 

We  repeated  the  scale  spectral  analysis  of  Part  2  for  other  data  sets,  the  97206-3  and  97206-15 
(reffered  to  as  03  and  15  in  the  captions).  The  results  are  quite  similar  and  agree  quite  well  with 
the  analysis  of  Part  1,  for  the  same  data  sets.  The  estimated  slope  and  log  intercept  processes 
obtained  by  segmentation  and  filtering,  and  no  averaging  over  adjacent  blocks,  is  a  bit  choppier 
that  what  we  find  in  Part  1  but  the  main  features  are  quite  comparable. 

As  we  have  already  mentioned  in  the  course  of  describing  the  analysis  of  the  data  in  Part  2,  the 
main  overall  consequence  of  the  careful  removal  (minimization)  of  segmentation  effects  is  that  the 
estimated  slope  process  shows  more  variability  that  appears  to  be  intrinsic. 
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